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In the mixed state of an extreme type-II d-wave superconductor and within a broad regime 
of weak magnetic fields (H c i <C H <C H C 2), the low energy Bogoliubov-deGennes quasiparticles 
can be effectively described as Dirac fermions moving in the field of singular scalar and vector 
potentials. Although the effective linearized Hamiltonian operator formally does not depend on 
the structure of vortex cores, a singular nature of the perturbation requires choosing a self-adjoint 
extension of the Hamiltonian by imposing additional boundary conditions at vortex locations. Each 
vortex is described by a single parameter 9 that effectively represents all effects arising from the 
physics beyond linearization. With the value of 9 properly fixed, the resulting density of states of 
Dirac Hamiltonian exhibits full invariance under arbitrary singular gauge transformations applied 
at vortex positions. We identify the self-adjoint extensions of the solutions found earlier, within 
the framework of the linearized Hamiltonian diagonalized by expansion in the plane wave basis, 
and analyze the relation between fully self-consistent formulation of the problem and the linearized 
model. In particular, we construct the low-field scaling form of the nodal quasiparticle spectra 
which incorporates the self-adjoint extension parameter 9 explicitly and generalizes the conventional 
Simon-Lee scaling. In a companion paper, we also present a detailed numerical study of the lattice 
d-wave superconductor model and examine its low energy, low magnetic field behavior with an eye 
on determining the proper self-adjoint extension(s) of the linearized continuum limit. In general, 
we find that the density of quasiparticle states always vanishes at the chemical potential, either 
linearly or by virtue of a finite gap. The low energy continuum limit is thus faithfully represented 
by Dirac-like fermions which are either truly massless, massless at the linearized level (mass ~ H) 
or massive (mass ~ vH), depending on the mutual commensuration of magnetic length and lattice 
spacing. 



PACS numbers: 74.25.Jb, 74.25. Qt 

I. INTRODUCTION 

In conventional s-wave superconductors the BCS 
quasiparticle excitation spectrum is gapped and, in pres- 
ence of a vortex defect, the quasiparticles form a dis- 
crete set of Caroli-deGennes-Matricon states 1 . These 
states are bound below the gapped continuum and are 
described by wavefunctions that decrease exponentially 
away from the vortex position. Once finite density of 
vortices is induced by a magnetic field, as is the case 
in the mixed state of type-II superconductors, these dis- 
crete levels broaden into extremely narrow quasiparticle 
bands whose effect on bulk properties is rather modest, 
as long as the magnetic field H remains well below the 
upper critical field H C 2- The situation is entirely differ- 
ent in high temperature cuprate superconductors: there 
the quasiparticle gap has a d x 2_ y 2 symmetry and van- 
ishes at four nodal points on the Fermi surface. The 
excitation spectrum is gapless, with linearly vanishing 
density of states at the Fermi level. As a consequence, 
there are no strictly bound states near a vortex defect in 
such a superconductor—. Instead, the quasiparticle wave- 
functions exhibit a power law decrease away from vortex 
position and the spectrum is continuous^. Once a finite 
magnetic field is turned on and vortex lattice appears 
in the mixed state, the low energy quasiparticles are ex- 
pected to form broad bands which dramatically influence 
thermodynamics, transport and vortex dynamics of high 



temperature superconductors. 

After an early semiclassical approximation due to 
Voloviki, who suggested that quasiparticles spectrum un- 
dergoes a Doppler shift E^ — > E^ — v(r) • k due to the 
vortex-induced superflow v(r) at a given point r, the 
quantum theory of quasiparticles in the presence of vor- 
tex lattice evolved along two somewhat separate lines. 
On one side, the self-consistent Bogoliubov-deGennes 
(BdG) equations for ci-wave superconductor were solved 
numerically in continuum, both for a single vorte^ and 
for a vortex lattice^. Additionally, the numerical solution 
was also obtained for BdG equations of a tight-binding 
lattice d-wave superconductor with boundary conditions 
corresponding to a periodic vortex array 6 . In all cases 
it was clearly demonstrated that a single vortex is inca- 
pable of truly localizing d-wave quasiparticles and broad 
low energy bands are evident in the quasiparticle excita- 
tion spectrum of a vortex lattice. 

Solving BdG equations in the inhomogeneous mixed 
state, particularly for a d-wave superconductor with its 
extended low energy quasiparticle states, is a daunting 
task. It is natural to inquire whether a simpler, analytic 
description might be devised which will allow one to ad- 
dress important and realistic physics questions including 
disordered vortex lattice, thermal and quantum fluctu- 
ating vortices, etc., where numerical solutions arc cither 
too forbidding or too opaque to be useful. An important 
step along this path is the linearized version of the BdG 
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Hamiltonian introduced by Simon and Lee£. Since the 
homogeneous system at zero field can be described by 
four species of massless two-component Dirac fermions 
residing at nodal points on the Fermi surface, it was ar- 
gued that in low magnetic fields H c \ <C H -C H C 2, when 
the separation between vortices given by magnetic length 
I becomes much larger than the size of the vortex cores 
- set by the BCS coherence length £ - the properties of 
low energy quasiparticles to the leading order in £/Z can 
still be described by a superposition of four independent 
Hamiltonians at each node£& 

By devising an eponymous gauge transformation, 
Franz and Tesanovio^ recast the nodal BdG Hamilto- 
nian as that of a Dirac particle moving in a zero average 
magnetic field and subject to effective long-range scalar 
and vector potentials whose point-like sources are located 
at vortex positions. This approach revealed that nodal 
quasiparticles couple to vortices through a combined ef- 
fect of two long-range terms: the semiclassical Dopplcr 
shift must be accompanied by a purely quantum mechani- 
cal Berry phase, which attaches a exp(±z7r) = (—1) phase 
factor to a wavefunction of a quasiparticle encircling an 
hc/2e vortej&iS. The FT transformation has been widely 
used to directly address, by combination of general sym- 
metry arguments2ii2*i2ii£ and explicit analytic and nu- 
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various 



merical calculations^ 
aspects of quantum mechanics of nodal quasiparticles 
in presence of vortices and will serve as the departure 
point for our discussion^. We should stress, however, 
that the results and conclusions of this paper are com- 
pletely general and remain unaffected if instead of the 
FT transformed BdG Hamiltonian one uses more famil- 
iar framework of the magnetic translations group applied 
to the original BdG equations. A natural question is 
why should not one simply use the magnetic translation 
group (MTG) and forgo the FT transformed problem al- 
together? The answer is twofold: obviously, in the FT 
transformed case one is dealing with representations of 
ordinary Bloch translation group that are far more con- 
venient than those of MTG. More importantly, the FT 
transformation is custom tailored for an efficient extrac- 
tion of the low-energy, long-distance sector of the original 
problem, a prime feature in numerous physical situations. 

While the problem of d-wave nodal quasiparticles 
interacting with vortex defects is of considerable the- 
oretical interest in itself, its rich phenomenology and 
relevance for numerous experiments on cuprates promote 
it to the forefront of modern condensed matter physics. 
In high temperature superconductors, the low field 
regime (H c i <C H <C H C 2) of the vortex state covers 
a large portion of the H — T phase diagram and a 
large body of experimental work thought to offer clues 
regarding microscopic origins of high temperature su- 
perconductivity is available both for analysis and future 
refinements. In particular, the high resolution STS/STM 
measurements of tunneling conductance in YBCO and 
BSCCO, which is proportional to the quasiparticle local 
density of states (LDOS), reveal peaks near vortices 



that appear at energies 5.5 meV in YBCO and 7 meV in 

BSCCO and at magnetic fields ~ 6 Tesla^ 3 .. The origin of 

the LDOS peaks is still a hotly debated subject and sev- 
eral -^^^ 32.24. 25. 26. 27.28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39.40 

describing various microscopic physics inside and around 
vortex cores, such as spin- and charge- density waves 
(SDW and CDW), Mott insulator, d-density wave and 
others, have been proposed to explain these features. 
This structure in the quasiparticle LDOS appears related 
to the experiments on neutron scattering 4 ^ which hint at 
the development of an antiferromagnetic order parame- 
ter localized in the vicinity of vortex cores. A mean-field 
factorization oft — J model on a tight-binding lattice and 
allowing for both e?-wave superconductor and AF orders 
does indeed, for a certain range of parameters in the self- 
consistent solution, describe a d-wave vortex state with 
an AF order developing inside vortex cores 2 ^* 2 ^*^ 3 ^^. 
The details of the calculations, however, depend strongly 
on the parameters of the model and the actual paring 
terms included in self-consistent approximation. In fact 
it is not even clear that the interior of vortices can be 
reliably described within the mean-field framework since 
the suppression of superconductivity inside vortices is 
likely to lead to a formation of a strongly correlated 
state characteristic of the pseudogap region in the phase 
diagram of cupratesSi 3 ^ .. 

The linearized continuum formulation of the bulk prop- 
erties of nodal quasiparticles has the distinct advantage 
of being essentially independent of such details pertain- 
ing to the structure near vortex cores. Nodal BdG quasi- 
particle behaves like a massless Dirac fermion and can 
be thought of as "critical"—: the Feynman path inte- 
gral trajectories of such particles lack any characteristic 
lcngthscale. This "criticality" implies a certain degree 
of "universality"^ 2 ^ 3 .: one expects that an effective low- 
energy theory of such Dirac fermions can be constructed, 
which incorporates all the effects of intra-vortex physics 
on the bulk properties of the BdG quasiparticles in terms 
of few simple parameters. A detailed derivation of such 
an effective theory, describing the influence of the physics 
in vicinity of vortex cores on properties of nodal BdG 
quasiparticles in the bulk, is the main goal of this paper. 

Following FT, we start by representing the motion of 
nodal BdG quasiparticles in a vortex lattice by mass- 
less two-dimensional Dirac fermions in combined vector 
and scalar potentials 8 . Both vector and scalar potentials 
are singular near vortices, and we will show that reg- 
ularization of this singular behavior is required, which 
itself depends on the short-scale features of the physics 
near vortex cores. The final effective theory is essentially 
the canonical FT Hamiltonian supplemented by a bound- 
ary condition at each vortex. The variety of possible 
physical scenarios describing different types of core and 
near-core behaviors is represented by a single dimension- 
less parameter 9 £ [0, tt) that characterizes the boundary 
condition and is associated with each vortex. 9 being di- 
mensionless, such boundary conditions can in principle 
modify the energy spectrum already at the leading or- 
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der in (kpl)^ 1 - to this order— the value of 9 uniquely 
determines the full extent of the influence that the struc- 
ture in or near vortex cores has on the bulk properties 
of nodal quasiparticles^. We generalize the conventional 
Simon-Lee scaling to include the explicit dependence on 
9 and numerically determine the scaling function. This 
generalized scaling form for the linearized quasiparticle 
energy levels is: 




where £ is a universal dimensionless function, n and k 
denote the band number and the Bloch momentum of 
the FT transformed states, respectively and Vp/v^ is the 
bare anisotropy of the d-wave nodes in a zero-field case. 

We find three prominent behaviors of a d-wave super- 
conductor as H — > (while still H ^> H c i) and determine 
the associated self-adjoint extensions of the linearized 
continuum problem: the first is a gapless behavior with 
the Dirac nodes intact - although renormalized by the 
field - but with their number doubled relative to the zero 
field result of four. This solution bears considerable re- 
semblance to the one reported earlier in Ref. 0. Second, 
we find new gapped spectra with no zero energy states. 
This self-adjoint extension produces a Dirac mass gap 
which scales as 1// and is thus a part of the leading or- 
der H — > scaling function, unlike the gaps of order l/l 2 , 
which vanish in the leading order scaling and are typically 
present even for the "gapless" case. The gapped solution 
seems to be related to the "interference" gapped spectra 
discussed in Refs. Il0ll47t Finally, we find 9 for which the 
quasiparticle bands lack (E, k) — * (— E, k) symmetry and 
for which the density of states is finite at the Fermi level - 
such solutions, however, at least for the superconducting 
gap parameter A not exceedingly small in comparison to 
hopping t, as is appropriate for near optimally doped or 
underdoped cuprates, are never found in our numerical 
calculations on the lattice ei-wave superconductor. These 
tight-binding calculations, reported separately^, always 
find either a linearly vanishing or gapped density of states 
at the Fermi level. Thus, the so-called "Volovik effect," 
a defining feature of an early semiclassical approach, is 
effectively absent from a fully quantum mechanical solu- 
tion. It would be important to verify this feature of our 
results experimentally, for example, by precision mea- 
surements of a specific heat or thermal conductivity. 

Alternatively, the above regularization and associated 
self-adjoint extensions can be implemented by introduc- 
ing a fictitious potential which suppresses the wavefunc- 
tions of nodal quasiparticles near vortices. This po- 
tential represents the nonlinear corrections to the lin- 
earized BdG Hamiltonian. Because of the well-known 
Klein paradox for Dirac particles, this potential cannot 
be chosen as a large scalar potential barrier, and instead 
must be realized as a position-dependent mass^, which 
grows rapidly as one approaches a vortex location but 
vanishes elsewhere. We demonstrate that the two de- 
scriptions of the regularized FT Hamiltonians are equiv- 



alent and yield essentially identical results in numerical 
calculations. The realization that the boundary condi- 
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FIG. 1: A magnetic unit cell containing two vortices. Panel 
(b) shows contours used in the application of the Stokes the- 
orem discussed in section ITTTl 

tions at the vortices must be externally imposed resolves 
the difficulty associated with previous attempts to find 
the spectrum of the quasiparticles numerically. Despite 
the gauge invariance of the linearized problem, the ex- 
plicit form of the FT Hamiltonian depends on the par- 
tition of the original vortex lattice into two sublattices 
A and B. Such choice is clearly not unique, and differ- 
ent selections are related by singular gauge transforma- 
tions. Obviously, all measurable physical quantities, ex- 
emplified by the density of states (DOS), must be invari- 
ant under such transformations. Surprisingly, numerical 
spectra of the linearized FT Hamiltonian exhibit small 
but persistent difference for distinct choices of A and B 
sublattices, when the wavefunctions are sought as linear 
combinations of plane waves^fi or when the Hamiltonian 
is solved by discretization on a real space mesbi!^. In 
what follows, we will demonstrate that the difficulties as- 
sociated with the linearized model indeed are not simply 
a nuisance related to the numerical procedure but rather 
reflect a genuine dependence of the bulk properties of 
quasiparticles on the internal structure of vortices. Even 
though the bulk of 2D CuO planes is dominated by pure 
nodal d-wave quasiparticles, the internal vortex structure 
affects the spectrum through boundary conditions that 
should be imposed at the vortex locations, due to sin- 
gular nature of the perturbation introduced by vortices, 
and the resulting energy spectra for different vortex core 
types are generally distinct to the leading order in Hvf/1- 
We reanalyze the solutions found earlier, and find that 
they correspond to mutually different sets of boundary 
conditions, "spontaneously" selected by a numerical pro- 
cedure used as well as the choice of A and B sublattices. 
Once the fixed boundary conditions are externally im- 
posed, the numerical procedure that incorporates them 
generates results that are fully gauge independent. 

The paper is organized as follows: in the next two sec- 
tions we introduce the linearized FT Hamiltonian and 
discuss the necessity of imposing boundary conditions 
(BC) at vortex locations. We show how these BC are 
intimately related to the self-adjoint extensions of the lin- 
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earized BdG HamiltonianAS. and discuss such extensions 
in detail in section 11111 After analyzing the symmetry 
properties of the linearized problem with and without 
BC in section IIVI we describe in section [V] a new pro- 
cedure used to address the problem, which incorporates 
the BC numerically. Section I VII establishes an alterna- 
tive framework for regularization of the linearized prob- 
lem using mass potentials and discusses the connection 
between the full, non-linearized BdG equations and self- 
adjoint extensions of the linearized model. This is fol- 
lowed by a general discussion of the relation between the 
d xy and d x 2_ y 2 continuum models in section lYlII Finally, 
we conclude with a brief summary of our results. 



II. THE LINEARIZED HAMILTONIAN 

It was argued^ that in the limit of weak magnetic 
field H < H c2 (but still H > H cl ), the BdG Hamilto- 
nian describing quasiparticles in an ideal vortex lattice 
can be formally expanded^ in powers of (fc^O -1 where 
magnetic length I is defined through the flux quantum 
$o = ^c/e as ( = y^o/H. In extreme type-II super- 
conductors, the regions separated from the vortices by 
distances sufficiently larger than vortex core size ~ £ are 
accurately described by the center-of-mass superconduct- 
ing order parameter A(r) = Aoe l< ^ r ' with a nearly uni- 
form amplitude Ao. One can then avoid solving the full 
self-consistent problem by eliminating the effective cou- 
pling constant in favor of Ao. In making this useful sim- 
plification, one implicitly excludes small disks of radius 
~ £ around each vortex in a 2D plane, where the am- 
plitude of the order parameter varies appreciably. This 
amounts to a tacit assumption that, in the limit £ <C I, 
the effects of the detailed structure around vortex cores 
on properties of the far-away nodal states in the bulk re- 
gions are entirely negligible. It is this assumption that 
needs to be revisited, as we will show. Initially, however, 
we will follow this assumption in order to understand the 
difficulties it generates and to introduce the notation. 

One starts from the BdG Hamiltonian for a d-wave 
extreme type-II superconductor: 
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The gap operator A is given by 

A = wp{p*> iPy> A W» " i^ A ( r )({ 5 - W)> ( 3 ) 

' F F 

where A(r) is the superconducting center-of-mass com- 
plex gap function, </> is its phase and curly brackets denote 
anticommutation operator {a, b} = (ab + ba)/2, accom- 
panied by the x y symmetrization. Note that the first 
termZ in Eq. © must be accompanied by the second^i^ 
to maintain gauge invariance of the BdG Hamiltonian - 
this amounts to replacing ordinary derivatives in the first 



term by the covariant ones: V — > V + (i/2)(V0), insur- 
ing minimal coupling of </> to the external electromagnetic 
gauge potential, with charge equal to 2e. Although the 
gap operator A in the above form has a d xy symmetry, 
the results for d x 2_ y 2 can be readily obtained after a sim- 
ple rotation by 7r/4. 

In a uniform external magnetic field H and in the pres- 
ence of a periodic array of superconducting vortices which 
the field induces in A(r), both diagonal terms and the gap 
operator in @ are invariant under magnetic translation 
group transformations (MTG) rather than being sim- 
ply periodic in space. Consequently, the eigenfunctions 
of Ti-BdG are not the ordinary Bloch functions but are 
instead the so-called magnetic Bloch functions and can 
be classified according to representations of MTG 50,5152 . 
An essential aspect of the problem is the observation that 
a simple Bravais lattice of superconducting vortices con- 
tains a magnetic flux equal to hc/2e per unit cell - this is 
just a reflection of the fact that a vortex is a topological 
defect in the phase of the gap function of Cooper pairs 
which carry charge 2e. In contrast, all representations 
of MTG in terms of single-valued eigenfunctions must 
contain an integer number of the electronic flux quanta, 
hc/e, per unit cell. It is therefore necessary to choose 
the unit cell of the MTG so it is at least twice as large 
as that of the vortex lattic o 53 i 54 . Apart from this condi- 
tion the shape and the size of such unit cell are arbitrary 
and different choices correspond to different "gauges" , i.e. 
different representations of MTG. The eigenfunctions as- 
sociated with these representations are different but the 
spectrum of eigenvalues of H.BdG, as measured by its den- 
sity of states, is the same for all these choices of MTG. 

FT observed that the above features of MTG can be 
used to recast the original problem in a more convenient 
form. One first partitions the original vortex lattice into 
two sublattices A, B (Fig.^) and then performs a singular 
gauge transformation Ti! = U~ 1 H.U: 
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e 



(r) = A (r) + ^(r) , (4) 



where 4>a(b) are contributions to the phase of the order 
parameter from vortices of A(B) sublattice (see^ and 
Appendix A). After the transformation the Hamiltonian 
assumes a periodic form: 
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where the transformed gap operator is given by 

A,, 



h 2 k 2 F 



{Px + a Xl p y + a y } (6) 



and two superfluid velocities are defined as 

mv A , B (r) =fiV^,fl(r)-eA/c 
a = (mvA — mve)/2 . 



(7) 

(8) 
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Properties of the superfLuid velocities and vg are 
discussed at length further in the text and in the Ap- 
pendix. The crucial feature is that their chirality, set 
by V x Vx(s)(r), vanishes on average and thus ~Va(b) 
are truly periodic functions in space. Consequently, the 
eigenstates of the FT transformed BdG Hamiltonian are 
ordinary Bloch waves that can be classified by their crys- 
tal momentum and band index. In the limit of weak fields 
it is reasonable to view the low energy wavefunctions as 
"perturbations" of the nodal states at H = 0, and there- 
fore the wavefunctions describing states close to the nodal 
point kf = (kp, 0) can be described as 

#(r) =e tkFX iP{r), 

where ipiijc) changes weakly on distances of order 
Under assumptions |V^(r)| <S fcp|^>(r)|, mv(r) <C hkp, 
and h\7 2 tp <§C kp\mv(r)\, the Hamiltonian to the leading 
order is& 

Uft = v F (Px + a x )<7 3 + VA<Ji{p y + Oy) + mv F v x , (9) 

where Oi are the standard Pauli matrices and v = (v^ + 
Vs)/2 is the gauge invariant superfluid velocity. The 
corrections to © fall into two categories: more familiar 
ones, of higher order in {kpl)~ x *C 1 become arbitrarily 
small as H — > 0. In addition, however, there are correc- 
tions governed by (kp^)^ 1 <C 1 even as H — > 0, which 
reflect the fact that nodes are located at finite momen- 
tum kp and describe rapid oscillations in quasiparticle 
spectra due to internodal interferenceiS^^I. 

The linearized Hamiltonian TLft describes a massless 
Dirac particle with anisotropic dispersion in presence of 
an internal gauge field a(r) and a scalar potential v x given 
by the component of v(r) along the node. This scalar 
potential is nothing but the Doppler shift of nodal quasi- 
particle energies in the superflow field produced by vor- 
tices. The internal gauge field a(r) - its minimal coupling 
to BdG quasiparticles betraying its topological origin - 
corresponds to an array of ±5 Aharonov-Bohm fluxes lo- 
cated at vortex positions and generates the required ±7r 
Berry phases for fermions circling around vortices&iS*^. 
This intrinsic topological frustration is entirely absent in 
semiclassical approximations™ but plays an essential role 
in the quantum mechanics of nodal quasiparticles in pres- 
ence of vortices. Note that a(r) is explicitly time-reversal 
invariant since the Aharonov-Bohm effect remains unaf- 
fected by the exchange of direction of a half-flux. This 
is an important feature and the necessary condition for 
invariance of physical results under different FT gauge 
transformations. Different FT transformations corre- 
spond to different choices of A and B sublattices (always, 
however, containing the same total number of vortices so 
that V x va(b) ( r ) remains zero on average) closely resem- 
bling different choices of the unit cell of the MTG. Under 
such transformations the Hamiltonian @ retains its gen- 
eral form but with half-fluxes in a(r) switching their signs 
from positive to negative depending on whether they be- 
long to A or B sublattice, while v(r) remains the same 
in all FT gauges. 



Hamiltonian Hpp is periodic and has a well defined 
Fourier representation. While the fact that both effec- 
tive potentials a(r) and v(r) diverge as 1/r near vortices 
is an obvious point of concern, resulting in matrix ele- 
ments of perturbation theory that decrease only as 1/k 
in momentum space, one still reasonably expects that 
the path toward finding the set of eigenfunctions of Hpp 
starts by simply expanding them in a plane wave ba- 
sis. The numerical solution obtained in this manner, 
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FIG. 2: Comparison of energy spectra for different choices of 
A and B sublattices. The spectra are similar at low energies 
but differ significantly as energy increases—. Lower left: the 
first magnetic Brillouin zone. 

shown in Fig. [5] was found to result in strongly dis- 
persive bands at low energies. The low-energy portion of 
the spectrum appears qualitatively similar to the zero- 
field case and, most importantly, preserves nodal point: 

e(k) = h^J v'pk 2 + v'£ky. The only significant effect of 

finite field is the renormalization of the original veloci- 
ties vp and va- The band structure calculation was also 
confirmed in£, where the spectrum and the eigenstates of 
Hpp Q were obtained by discretization of the Hamilto- 
nian in real space. The persistence of nodal points was 
further demonstrated-*^ to be a consequence of the gen- 
eral symmetry properties of Hpp. 

As noted previousljii-, the choice of A-B sublattice ar- 
rangement is not unique in the FT transformation, being 
merely a choice of singular gauge. Clearly, measurable 
quantities such as the density of states (DOS) should not 
depend on the partitioning of vortex lattice into various 
A and B sublattices. On the other hand, the Hamil- 
tonian Hpp has a formal dependence on the FT singu- 
lar gauge choice, and surprisingly, the numerical solution 
exhibits small but persistent differences for different ar- 
rangements of A and B sublattices (Fig. |2J). The dis- 
crepancy remains regardless of whether the computation 
is done by expanding the solution in plane wave basis&iS 
or by discretization in real space-. Note that these dif- 
ficulties cannot simply be attributed to the singular na- 
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ture of the FT gauge transformation: the energy spec- 
trum of the Hamiltonian in the original representation, 
where the eigenfunctions are chosen as magnetic trans- 
lation group eigenstates, has DOS which also exhibits 
a similar discrepancy for different choices of a magnetic 
unit cell (e.g. rectangular vs. square). As a remedy, the 
tight-binding formulation of a <i-wave BdG Hamiltonian 
and FT transformation was introduced inAfi, where it was 
verified that this so-called "ABAB vs. AABB" problem 
does not appear and the spectra are explicitly invariant 
under singular gauge transformations for arbitrary size 
of the tight-binding mesh. In a way, this is the resolution 
of the whole "ABAB vs. AABB" problem: Cu0 2 planes 
in high temperature superconductors are actually well 
described by an effective tight-binding Hamiltonian (of a 
t—J, Hubbard or a related variety) and the lattice d-wave 
superconductor and its BdG equations is precisely what 
we should be considering. The problem with the sin- 
gular 1/r behavior never arises since BdG quasiparticles 
live on sites while vortices are "located" in the interior 
of plaquettes of a tight-binding lattice. 

Nevertheless, the appeal of the continuum formulation 
is undeniable. By representing nodal quasiparticles as 
massless Dirac fermions one is hoping that their inter- 
actions with vortices can be described by a relatively 
simple and elegant effective continuum theory which can 
then become the starting point for analytic exploration 
of more complex problems involving fluctuating or dis- 
ordered vortices, thermal and quantum fluctuations of 
vortex-antivortex pairs, etc. It is therefore highly de- 
sirable to understand the source of difficulties besetting 
the continuum description of the linearized model and to 
overcome them£&. Moreover, as pointed in the Introduc- 
tion, a detailed structure of vortex cores in cuprates is 
unknown at present and all possible tight-binding ap- 
proaches will necessarily suffer from being dependent 
on short-range physics details, which, as one suspects, 
should have only a limited effect on the nodal quasipar- 
ticle wavefunctions in the bulk and thus on any effective 
low-energy theory. The linearized theory, in contrast, 
is formulated in terms of a universal Hamiltonian J^}, 
which is completely independent of short-range details, 
including the structure of vortex cores. 



III. SELF-ADJOINT EXTENSIONS 

The origin of the difficulties with different choices of 
FT singular gauge can be traced back to the assump- 
tions made during the linearization procedure outlined 
above Eq. ©. Although the conditions necessary for 
linearization are well satisfied for weak magnetic fields 
and in the bulk, far away from vortex cores (r ~ I), 
they are violated at distances r ~ £, since the neglected 
terms diverge as 1/r 2 near vortices. Thus one suspects 
that additional regularization might be required to take 
into account the effect of such terms. Let us nevertheless 
assume for a moment that the dependence on A-B ar- 



rangement is caused by a numerical procedure. Since the 
scalar and vector potentials a(r) and v(r) diverge as 1/r 
as one approaches a vortex, one naturally expects that 
by finding the asymptotic behavior of the eigenfunction 
near vortices, one can use it to his/her advantage to cure 
the singularities. The asymptotics in fact were found an- 
alytically by Melniko\*£&. For brevity, in this article we 
will consider mostly the isotropic case (vf = «a)| the 
results are easily generalized to the anisotropic case. 

Close to a vortex two linearly independent solutions 
of the eigenvalue problem for the linearized Hamiltonian 
TLft are described by the wavefunctions that diverge 
near the vortex as r -1 / 2 . The angular dependence of the 
divergent part of the wavefunctions y/rtf>(RA + r) near 
vortex A is described by 

(he * ^ iv _ 1 j+C 2 e ■ [ ie -i v + 1 ) , (10) 

where <p (not to be confused with 0(r), the phase of the 
gap function) is the polar angle around a particular vor- 
tex. The above is only the dominant part of each eigen- 
function as r — > 0; the prefactors of this diverging term 
as well as subdominant terms are different for different 
eigenfunctions. Generally, the weight of the divergent 
term in an eigenfunction decreases as its energy eigen- 
value increases. The angular dependence of the singular- 
part of the wavefunctions near vortices of B sublattice 
y/rip(R,B + r) is given by an expression similar to l(TT)jl : 

C,e > [i_ leiv )+C 2 e . [ i + le iv)- (11) 

The asymptotics above can be checked by retaining 
only the divergent parts of the potentials v and a in the 
FT equations with the radial dependence proportional 
to 1/r. Such "zeroth-order" approximation describes a 
single vortex, since both the superfluid velocity and the 
vector potential can be represented as a sum (see Eq. 
(|A9I) of the Appendix). Each term of the sum represents 
contribution of an individual vortex, and by neglecting 
all terms of the sum that are finite only the contribu- 
tion of a single vortex is taken into account. Then, one 
can easily verify that (|l(Jfl or (|11|) . for a vortex belong- 
ing to class A or B respectively, are exact eigenfunctions 
corresponding to zero energy. Returning now to the lat- 
tice problem, the non-singular part of the Hamiltonian, 
V reg , which is due to other vortices, can be treated as a 
perturbation. Such finite terms do not modify the lim- 
iting small r behavior of the zero energy eigenstates (if 
there are any). Moreover, under the action of perturba- 
tion V regi the nonzero energy eigenstates are mixed with 
the singular wavefunctions IjlOfl or (|llfl . and therefore 
are described by the same limiting behavior close to any 
particular vortex. 

The eigenfunctions of H.ft obtained numerically by 
the plane wave expansion close to a vortex will be shown 
below (see Fig. [3J to have asymptotic behavior given 
by either (|10|) or l|ll|) for vortices belonging to sublat- 
tice A or B, respectively. Additionally, the increase of 
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the wavefunctions according to power law r -1 / 2 and an- 
gular dependence given above near vortices at distances 
£ <C r <C I has been explicitly verified 60 for the solutions 
of BdG equations describing the tight-binding lattice d- 
wave superconductor^^. 

Note that Eq. (|1U[) describes two linearly independent 
solutions both of which diverge as 1 /y/r as one approaches 
the vortex. For non-singular potentials only one of the 
solutions would have been square-integrable, and the sec- 
ond solution would have to be discarded. Thus, before 
one proceeds, one must decide whether such divergent 
solutions are permitted in the spectrum. 



A. Dirac particles in the field of Aharonov-Bohm 

flux 



potentials and a more careful analysis is needed. The do- 
main of finite wavefunctions should be enlarged to include 
a specific linear combination of l/y^-divergent solutions. 
Using the language of functional analysis, self-adjoint 
extensions of the Hamiltonian should be constructed 6 ^. 
Gerbert and Jackiw found that the mathematically al- 
lowed linear combination of divergent wavefunctions is 
not unique but forms a one-parameter family 

x oc smw y ieiVJi/2 ( Er ) j + co&y y_ ie ^j_ i/2 ( Er) J 

— ( s[li9 \ 



This issue is not entirely new, and similar problems 
were encountered in the context of cosmic stringaSl. 
There, a Dirac particle is considered in the field of a sin- 
gle Dirac string carrying a fractional magnetic flux. We 
will restrict our discussion to the relevant case of a single 
half-integer flux: 



Px 

Py 



a„ 



Py 

-Px 



= E 



(12) 



where a = is the vector potential corresponding to 

the Aharonov-Bohm half-flux. The Hamiltonian can be 
obtained from where the last term describing scalar 
potential is set to zero. After rotation xi = u — iv. Y2 = 
u + iv we obtain Hamiltonian in the basis used in 61: 



Hgj 



Px + a x - i(p y + CLy) 

Px + a x + i(p y + a y ) 



(13) 

The eigenstates of Hqj can be expressed through Bessel 
functions as 



ie lv J r 



1/2 

+3/2 



Co 



J_ 



-1/2 



3/2 



Gerbert and Jackiw noticed that for all angular chan- 
nels except n = — 1 the square integrability requirement 
specifies which of the two independent solutions in 114|) 
should be present. For n = — 1, however, there is an am- 
biguity as both solutions are square integrable, but diver- 
gent as l/y/r at the origin. They found that keeping both 
wavefunctions in the spectrum results in an overcomplete 
basis, with such pathologies as imaginary eigenvalues of 
the Hamiltonian. On the other hand, the requirement of 
continuity (non-divergence) of the wavefunctions at the 
origin is too strong and causes the loss of completeness 
in the angular momentum channel. 

The complications are due to the singular nature of the 
vector potential a(r) at the origin. The usual spectral 
theorems, normally derived for bounded symmetric (Her- 
mitian) operators, do not hold automatically for singular 



where 9 S [0, it) is a real parameter that labels distinct 
self- adjoint extensions. The parameter 9 expressing the 
boundary condition cannot be found from the model that 
treats the string as a "black box" ; it depends on the 
short-scale structure of the string. 

A string described by a divergent vector potential at 
r = is of course only an idealization. To find appropri- 
ate 9 one has to consider the physical regularization of 
the problem. The simplest case of magnetic field concen- 
trated in a thin cylindrical shell of small, but finite radius 
e when e — ► was considered in 1(341 By matching solu- 
tions inside and outside the core, the authors found that 
for radially extended symmetrical distribution of mag- 
netic field inside the core one obtains 9 = 0, implying 
that the lower component of spinor y stays regular at the 
origin. The procedure can be repeated for physical, ex- 
tended fluxes carrying arbitrary half-integer flux <!>, and 
it was found that for $ = 1/2,3/2,5/2,7/2. . . the self- 
adjoint extension is the same 9 = and the solutions for 
positive half-integer fluxes are related to each other by 
singular gauge transformations x( r ) ~ * exp(iiV<y9)x(r). 
Quite surprisingly, the negative extended half-integer 
fluxes $ = —1/2,-3/2,... belong to a different class 
9 = it/2 which corresponds to wavefunctions with regular 
upper components of x- Thus, the symmetry $ — ► $ + 1 
is broken around $ = 0. Due to the unbounded increase 
of the wavefunctions, a relativistic fermion can "probe" 
the short-scale structure of the core and the direction 
of the flux is "exported" to the exterior through bound- 
ary condition 64 . Although physical, extended fluxes do 
have different limiting solutions for opposite orientations 
of magnetic field as the size of the core is taken to zero, 
one can still perform singular gauge transformations, of 
course, provided that the boundary condition 9 is cho- 
sen correctly to reflect the direction of the real, physical 
magnetic field. 

A necessity of considering self-adjoint extensions is not 
restricted to a single-flux problem. In the case of two 
fluxes the problem can also be solved essentially exactly; 
the solution is presented in the Appendix B. 
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B. Dirac particle in the field of periodic array of 
singular scatterers 

The above analysis of a single flux problem^ is facili- 
tated due to the knowledge of exact eigenfunctions allow- 
ing for application of the von Neumann theorem on self- 
adjoint extensions^. To consider an array of vortices, 
we will rederive the result l|15[) by following a different 
procedure, which will allow us to find the self-adjoint ex- 
tensions corresponding to the vortex lattice problem, for 
which exact eigenfunctions are unknown. We consider a 
general Hamiltonian 



H 



Px + Vll Py + Vl2 

Py + V 2 l ~Px + V 2 2 



(16) 



where Vij are arbitrary periodic functions which can di- 
verge at most as 1/r at a certain finite set of points in 
unit cell, such as positions of fluxes or vortices. 

In order for Hamiltonian J5J to be a symmetric opera- 
tor (i/j*Hx) = (x*Hil>)*, the condition 



[dxiXi^i ~ X2^ 2 *) + dyixi^ + X2^i)] d 2 r = (17) 



must be fulfilled. This is satisfied automatically if x( r ) 
and ip(r) have different crystal momenta, and therefore 
we concentrate on the subspace of the wavefunctions with 
the same crystal momentum k, but different band indices. 
Symmetric operator is self-adjoint if the domain D{H) 
coincides with that of its adjoint D(H*). In other words, 
we have to find such a boundary condition on ?/;(r) so 
that the adjoint wavefunctions %(r) satisfy precisely the 
same condition. For the subspace of functions charac- 
terized by the same crystal momentum k, the integrand 
of l|17fl is periodic in space, and therefore the integral 
over entire space can be replaced by an integral over one 
unit cell. It may appear that the integral automatically 
vanishes; by Stokes theorem it can be transformed into 
a contour integral along the edges of the unit cell, which 
in turn vanishes due to the periodicity of the integrand. 
Indeed, this conclusion is valid for regular wavefunctions. 
It must be remembered, however, that the necessary con- 
dition for Stokes theorem to hold is the continuity of the 
partial derivatives in the integrand in (|17fl inside the con- 
tour, and in our case, with ip tx l/y^ near a vortex, the 
theorem does not automatically apply. To proceed with 
the analysis, we define 



A y = X21P2 - X11P1 



(18) 



and separate the contribution of singularities near the 
vortices by drawing discs I\ of radii e$ surrounding each 
vortex (see Fig.^b) for the case of two vortices per unit 
cell). 

After Stokes theorem is applied to the exterior of the 
discs within the unit cell, where the wavefunctions are 
finite, the contribution from the edges cancels by period- 
icity, and we find that the sum of the contour integrals 




FIG. 3: Typical eigenfunctions of the linearized equation 
found by expanding the wavefunctions in plane wave basis 
without explicitly imposing boundary conditions at vortex lo- 
cations. N = (2n+ l) 2 plane waves are included in numerical 
solution. Right panels: y'r|u(r)| and xA'I'H 7 *)! along fixed di- 
rection ip = as functions of distance from the vortex. Each 
plot approaches a straight line as the accuracy is increased. 
The convergence is non-uniform exhibiting Gibbs overshoot 
as expected, because of the singularity of the wavefunction at 
the origin. Upper left panel: ratio \u(r)/v(r)\ used to extract 
6. The intercept of the limiting linear dependence (0.3+0.67r) 
is used to find 6 as shown in the lower left panel. 



along the boundaries of the discs Ti and the area integrals 
over the interior of Ti equals 

V" / (sin tp A y (ej, tp) - cos(pA x (ei, (p))dip = , (19) 

where i labels the vortices in unit cell and ip is the polar 
angle around the ith vortex. Obviously, if one demands 
that 7/j(r) are regular at the origin, then the adjoint do- 
main of x includes all solutions that behave as 1/ \pr , 
leading to D{H) ^ D(H*). Now note that in the limit 
€i — > only the singular part of the wavefunctions (|1C)|) 
and (|11|) contributes. Since the limits €i — > can be taken 
independently, it is straightforward to find from <|19f) that 
the self-adjointness requirement D(H) = D(H*) fixes 
the relative phase between two asymptotic solutions 62 . 
When applied to the array of fluxes we find that this re- 
uircmcnt leads to the same boundary conditions as in 
1. In terms of (u,v) these boundary conditions in the 
vicinity of fluxes are 



-i<j> 



$0 

for +— fluxes 
(20) 

$0 

for fluxes . 

2 

(21) 



9 



These boundary conditions are actually equivalent to 
those of Ref. El 

For the present problem of d-wave quasiparticles in a 
vortex lattice the above requirement of self-adjointncss 
D(H) = D(H*) translates into the following boundary 
conditions near A vortices: 

s/rip(R A + r) -> 

cosOe * { ie _ iv _ 1 )+ S mOe * [ ie -i v + 1 )' 

(22) 

The boundary conditions at B vortices differ by the over- 
all phase factor exp(z</?): 

s/fip(R B +r) -> 
cosOe * ^._ eitp )+sm9e * ^ + ^ j. 

(23) 

Thus, if the interior of a vortex is treated as a black 
box, in order to fully specify the problem, the linearized 
Hamiltonian Hpx must be supplemented by boundary 
conditions at the locations of vortices^. Every boundary 
condition depends on a single parameter 9, which is easily 
shown to be independent of an A-B assignment. Indeed, 
the transformation from one choice of A and B sublattices 
to another is given by a unitary matrix 

/ e i<t>A-i<i>A' n \ j. j. 

II a a, — I , I — P l( t>A-i<l> Al . -i 

Uaa ~ I o e ^B-^B J ~ e 1 ' 

since (ft A + 4>b — 4>A' + <\>B' — 4>{ v )- Uaa 1 is proportional 
to unit matrix, and the asymptotic behavior of the wave- 
functions, specified by 0, is clearly invariant under all 
such transformations. Note that parameter 9 for a given 
vortex has been intentionally defined in this particular 
way for convenience. It is independent on whether this 
individual vortex belongs to an A or a B sublattice, and 
is therefore a scalar under FT singular gauge transfor- 
mations. Naturally, one must recall here that the actual 
form of the boundary condition does change according to 
Eqs. H22I23|) depending on whether a particular vortex 
belongs to an A or a B sublattice. Furthermore, note 
that there is no a priori requirement that parameters 9 
are equal for different vortices within unit cell. Just as 
in the case of a single Dirac string, parameters 9 cannot 
be determined from the linearized Hamiltonian. Rather, 
the boundary conditions 9 are determined by the short- 
ranged physics of fully self-consistent BdG equations and 
include all effects of higher orders in (kpl) -1 expansion. 

Equipped with this new understanding we now revisit 
the "ABAB vs. AABB" problem and consider the numer- 
ical computations performed earlier^iifi. There boundary 
conditions were not explicitly enforced and the numerical 
procedure itself "spontaneously" selected a particular set 
of boundary conditions. In calculations done inlMlol the 
wavefunctions are represented as linear combinations of 



plane waves, and, after a substitution into the linearized 
FT Hamiltonian, the series is truncated at some number 
of terms N = (2n + I) 2 , which is gradually increased until 
convergence is achieved: 




Rather than analyzing the spectra we first concentrate 
on the wavefunctions. The solutions ip(r) = (u(r),v(r)) 
are shown in Fig. [3] The wavefunctions indeed approach 
f/v / r dependence on the distance from the vortex as 
the number of reciprocal lattice vectors included in the 
solution is increased. The values of 9 for ^?„(k) can 
be extracted from the wavefunctions by various meth- 
ods. As an example, the procedure involving ratios 
linip^o Wip,^ — 0)\/\v(p,<p — 0)| is illustrated in Fig.0 
Since there are two vortices per unit cel]»> we label them 
"A" and "B" to indicate their belonging to two different 
sublattices - to avoid any confusion, it is important to 
stress here that the labels "A" and "B" refer to two phys- 
ically distinct vortices rather than to an FT gauge label 
attached to an individual vortex (the quotation marks are 
used to highlight this difference). The value 9 for vortex 
"A" turns out to be close to 9 = (0.75 ± 0.03)7r, while 
vortex "B" has asymptotic form with 9 = (0.25 ± 0.03)7r. 
The conclusion remains valid not only for the vicinity of 
nodal points, but also for all bands at arbitrary k for 
which 9's could be extracted. As discussed in a separate 
paper these values of 9's are among several consistent 
with the general symmetry requirements that must be 
obeyed by the solutions of the full BdG equations. They 
thus faithfully represent a possible solution for the low 
energy spectrum of BdG equations. The meaningful de- 
termination of 9 for high energy bands was beyond the 
precision of our calculations since the contribution of the 
singular part of the wavefunctions decreases for higher 
bands, and one has to consider distances that are ex- 
tremely close to the vortex to extract the asymptotic 
behavior. This, in turn requires diagonalization 
of matrices with unrealistically large N. 

We then performed a similar analysis of several other 
arrangements of A and B sublattices. We found that 
the eigenstates of "AABB" and "ABAB" lattices for 
four vortices per unit cell and similar vortex lattices ro- 
tated by 7r/4 with respect to the nodal directions also 
exhibit the same pattern: the plane wave expansion pro- 
cedure drives "A" vortices to ~ 0.757T and "B" vor- 
tices to 9 ss 0.25-7T. The discrepancies encountered earlier 
that gave rise to the "ABAB vs. AABB" problemiS are 
now easy to understand, since the problems solved with- 
out externally imposing boundary conditions at vortex 
sites corresponded to physically different situations once 
boundary conditions were "spontaneously" generated by 
the chosen procedure, i.e. the two problems solved were 
not connected by an FT singular gauge transformation. 
In particular, the "ABAB" sublattice partition corre- 
sponds to the boundary conditions at vortex locations 
with 9 = ±7r/4 assigned in a checkerboard arrangement. 
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In contrast, if vortices "A" and "B" are chosen to form 
series of parallel lines, as in the "AABB" partitioning, 
the numerical procedure spontaneously selects boundary 
conditions arranged as interchanging lines of 9 = tt/4 
and 9 — — it/ 4. Obviously, these are two genuinely dif- 
ferent physical situations in light of the condition that 9 
for a particular vortex must remain unchanged under FT 
gauge transformations. 

It is clear that the boundary conditions at vortex lo- 
cations must be fixed separately and changed appropri- 
ately under FT transformation following Eqs. 122123( 1 . as 
anticipated inllOl A specific 9 assignment can be imple- 
mented by a procedure which represents a modification 
of a plane wave expansion similar to the orthogonalized 
plane waves method. This procedure is introduced and 
discussed in detail in section V and yields results in mu- 
tual agreement for different FT singular gauge choices. 
Alternatively, the boundary conditions can be emulated 
by adding a mass term to the linearized Hamiltonian that 
grows rapidly as one approaches a vortex core. The latter 
approach is described in section VI. 



Group element h 


Transformation of v and a 


Mirror symmetry m x : 
m x r = (1/2 - x,y) 


v x (gr)=v x (r) v y (gr) = -v y (r) 
a x (gr) = a x (r) a y (gr) = -a y (r) 


Mirror symmetry m y : 
niyV = (x,l/2 - y) 


Vx{gr) = -V x (r) Vy(gr) = v y (r) 
a x (gr) = -a x (r) a y (gr) = a y (r) 


Inversion / = m x m y : 
Ir = 2Ra - r 


v(gr) = -v(r) 
a(gr) = -a(r) 


Inversion A — > B 
Pv = -r 


v(gr) = -v(r) 
a(gr) = a(r) 


Pm x r = (x + 1/2, -y) 


v x (gr) = -v x (r) v y (gr) = v y (r) 
a x (gr) = a x (r) a y (gr) = -a y (r) 


Pm y r = (-x,y + l/2) 


Vx(gr) = v x (r) vy(gr) = -v y (r) 
a x (gr) = -a x (r) a y (gr) = a y (r) 


Translation s — PP. 
sr = r + 2Ra 


v(gr) = v(r) 
a(#r) = -a(r) 



TABLE I: Symmetry properties of superfluid velocity v and 
internal gauge field a. 



IV. SYMMETRIES OF THE SINGLE-NODE 
HAMILTONIAN 

In the previous section we showed that the linearized 
FT Hamiltonian Hft must be supplemented by bound- 
ary conditions at vortex locations. These boundary con- 
ditions are fixed by a dimcnsionless parameter 9 € [0, it). 
The specific value of 9 that should be used is determined 
by the physics unfolding at the lengthscales shorter than 
those included in the linearized description. Thus, we are 
at an impasse; we apparently must solve for the physics 
beyond linearization to fully specify the linearized prob- 
lem itself. We take the initial step toward such a solu- 



tion in a separate paperjS 7 . where we study a lattice d- 
wave superconductor. In the present paper, we continue 
our analysis of different 6*'s by focusing on the symme- 
try properties of Hft hi combination with those of the 
original, non- linearized BdG problem I|2I5|I . 

Certain aspects of the symmetry properties of the 
linearized Hamiltonian Hpp were already considered 
previouslySii^. For example, it has been shown that the 
Dirac node of the zero-field problem is not destroyed by 
the inversion-symmetric vortex lattice; the only situa- 
tion considered, however, is the one with no boundary 
conditions imposed on the vortices. Here, we extend and 
generalize these early results by analyzing the complete 
group of symmetry transformations of Hft and examine 
the consequences of imposing the boundary conditions 
(|22I23(1 . For convenience, only the simplest choice of a 
unit cell (Fig.^) containing two vortices per unit cell will 
be studied. Using the Fourier representation l|A4(l and 
(|A5I) . the superfluid velocity v(r) and the internal gauge 
field a(r) can be shown to transform under geometric 
point transformations h according to Tabled Note that 
under inversion / around a vortex both v(r) and a(r) 
change sign, while the inversion operation P around the 
midpoint between vortex A and vortex B which inter- 
changes vortices A and B, leaves the internal gauge field 
a(r) invariant: a(— r) = a(r). 

The operations h do not form a group by themselves, 
since their product can generate a pure translation by a 
lattice vector. Rather, the group is formed by elements 
Tr/i where Tr, is an arbitrary translation by a lattice 
vector R. One can easily check that the group can be 
generated by forming products of primitive translations, 
reflections m x , m y and an inversion P. The symmetry 
of the potentials v(r) and a(r) is higher than the group 
just described, including also operations generated by ro- 
tations by 7r/4 around a vortex; these additional opera- 
tions, however, do not correspond to any symmetry of the 
linearized Hamiltonian due to the last term in @ reflect- 
ing the fact that the boost in the direction of one of the 
nodes breaks the four-fold symmetry of the original BdG 
Hamiltonian even in the isotropic case vf — ka- Only 
after the contributions of all four nodes are combined, 
will this symmetry be fully restored. 

It should be emphasized that the symmetry properties 
discussed in this section refer to the linearized Hamilto- 
nian Hft and its eigenfunctions ip(r). ^p(r) are related 
to the wavefunctions ^(r) in the original BdG basis ac- 
cording to 

(r) = e lkFX " V(r). 

w V e - l< M r ) J w 

All transformations introduced above are understood as 
infra-nodal. For example, inversion P transforms ip(r) to 
ip(— r), and the transformed wavefunction in the original 
basis is given by 

P*(r)=e^^ e q e _£ w U-') ■ (25) 
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Properties of the original non-linearized Hamiltonian (JSJ 
under various symmetry operations and their relation to 
the symmetries of Hft will be discussed separately in 
section IVT1 

Before analyzing the problem with imposed bound- 
ary conditions, we first consider formally the linearized 
Hamiltonian Hft with no boundary conditions enforced 
at vortex locations. Although such "unrestricted" prob- 
lem is in principle not guaranteed to represent any partic- 
ular physically meaningful linearization of the full BdG 
Hamiltonian (J2J), this discussion will serve as the basis 
for subsequent analysis. After the boundary conditions 
are enforced, the allowed symmetry operations will form 
a subset of the full transformation group of the "unre- 
stricted" Hamiltonian, which depends on the choice of 
the self-adjoint extension. 

Using the properties of superfluid velocities from Table 
[Jit is straightforward to verify that if is an eigenfunc- 
tion of Hamiltonian TIft with crystal momentum k and 
eigenvalue E then VE^, = D(g)^fk is an eigenfunction of 
H ft with crystal momentum and eigenvalue ±E accord- 
ing to the Table II I II in the Appendix. In identifying k' 
after the transformations involving phase factor 

/(r) = exp(-iS<p) = exp(i^ s - i<fi A ) 

it is important to keep in mind that /(r) is not peri- 
odic but anti-periodic - it changes sign under primitive 
translations along x and y directions (see Appendix): 

f(x + l,y) = -f(x,y) (26) 
f(x + l,y) = -f(x,y) (27) 

Transformations D(g) form ray representation^ of the 
group, and the standard methods can be applied to an- 
alyze the degeneracy of the eigenstatcs. In addition to 
geometric transformations and pure translations D(Tr), 
there is an additional anti-linear complex-conjugation 
transformation: 

L>(C)* = ie'^a^ 

D(C) commutes with all other operators up to a phase 
factor, and combined with the operations listed in the left 
portion of Table II I II from Appendix, C generates 8 new 
eigenfunctions with energy ±E. Since Bloch basis has 
already been chosen to diagonalize ordinary translation 
operators we have the total of 16 non-trivial transforma- 
tions, which generate 8 eigenfunctions with energy E and 
8 eigenfunctions with energy — E. 

If k is an arbitrary point of the Brillouin zone that 
does not lie on one of the symmetry lines, all 8 states 
with energy E correspond to distinct crystal momenta 
k' forming a "star" shown in Fig. 01 and are therefore 
linearly independent. Consequently, for general k the 
symmetry of the Hamiltonian is not sufficient to result 
in degeneracy of the bands at the same k. For every state 
(E, k) there is a state with the opposite energy (— -E^k) 
obtained by applying transformation CP to the origi- 
nal wavefunction and seven other pairs of states with 



energies ±E with crystal momenta k' belonging to the 
"star" . For certain particularly symmetric values of k 
two or more points of the "star" coincide, and there- 
fore the spectrum may become degenerate, as will be 
discussed shortly in detail after the boundary conditions 
are chosen. Note, however, that the Hamiltonian appar- 
ently possesses higher symmetry than what is reflected 
in the band structure spectrum shown in Fig. [5] which 
was obtained numerically by the plane wave expansion 
(without imposing boundary conditions explicitly near 
vortices In particular, the symmetry analysis in- 
dicates that the node at k = implies also the identical 
node centered around the corner of the Brillouin zone 
7r, while no such additional nodes are seen in numerical 
computation. In view of our previous discussion of the 
role of boundary conditions, the discrepancy is easy to 
understand: we have shown that, even if the boundary 
conditions at the vortices are not explicitly enforced from 
the outset, the numerical procedure selects a certain self- 
adjoint extension, which has the effect of reducing the 
symmetry of the problem. 

To study the effect of boundary conditions on the sym- 
metry of the spectrum, we note that, if one starts with 
Hamiltonian Ti. ft with boundary conditions on the wave- 
functions given by parameters #1(2) near the two 
vortices, then after general transformation wavefunctions 
D(g)ifj(r) have asymptotics described by a set of differ- 
ent boundary conditions #^ 2 )- ^he new parameters $i( 2 ) 
are shown in the rightmost column of Table II I II from 
Appendix. Note that the spectrum is symmetric under 
(k, E) — ► (k, — E) and only if 9\ = —8 2 - In the case 
when 6*i = B\ that leads to only two possible choices: 
01 = &2 = an d Oi = 02 = it/2. For these extensions, 
however, the inversion property E(k) = E(— k) is lost 
and we find that no self-adjoint extension exists with the 
spectrum resembling that of the zero- field Dirac particles, 
if all vortices are assigned the same value of 9. 

Plane wave extension. The numerical analysis of sec- 
tion ^H] showed that the straightforward plane wave ex- 
pansion "spontaneously" generates solutions with the 
boundary conditions #1(2) ~ ±7r/4. We now first demon- 
strate that the spectrum of the self-adjoint extension with 
&i(2) — ±7r/4 indeed has the symmetry properties consis- 
tent with the numerical solution and then clarify why the 
plane wave expansion procedure favors these particular 
boundary conditions. 

The subset of symmetry operations that leaves #1(2) = 
7r/4 invariant and its star of the k- vectors is shown in 
Fig. [SJ The special points in BZ where the spectrum 
can potentially become degenerate are T, X, Y, and M, 
where all four vectors of the star coincide^. Whether this 
indeed leads to degenerate eigenfunctions depends on k 
and should be analyzed in each case separately. Using 
Table IIIII from Appendix, the representations matrices 
D(PMy) = g x , D(CPI) = 32, and D{Cm x ) = g 3 satisfy 
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FIG. 4: First Brillouin zone is shown with the star of 
wavefunctions generated by applying symmetry operations to 
*&k(r). For convenience both, positive energy states (open cir- 
cles) and negative energy eigenfunctions (crosses) are shown. 



g\ = e ikyl 9\92 = ier lkyl gz gig?, = -ig 2 

gigl = -ie liK ~ k y )l g 3 g 2 gi = -e i{ - k * +k ^ 1 gi~gl,=ig\ 
g 3 gj = ie lkJ g 2 g 3 g2 = ie lk * l gi g 3 g^ = 1 

(28) 

At r and M-points g 2 g 2 = — 1 and consequently there 
are no one-dimensional representations. The lowest di- 
mension of irreducible representation of relations Q28JI. 
realized by Pauli matrices gi — <7j, is two. Thus, the 
eigenstates in the center and in the corner of the BZ are 
at least doubly degenerate. 

At points X and Y of Brillouin zone 

gl^(v) = *(r), 

and there are two one-dimensional representations g\ — 
±1, g 2 = ±iA, and g 3 = A, where A is an arbitrary 
complex parameter with unit norm |A| = 1. The energy 
bands therefore are non-degenerate at X and Y. These 
results are in complete agreement with the numerical so- 
lution obtained by Franz and TesanovicS which is shown 
in Fig. |U 

To understand why the boundary conditions with 
Q\(2) — ±7r/4 are favored by the plane wave expansion 
procedure, we recall that in this method the eigenfunc- 
tions of Hft are sought as 

Mr)= E f"^e l(Q+k)r (29) 

QSEjv \ V qJ 

where the set of reciprocal lattice vectors £ at in the ex- 
pansion form a square grid (2 A + 1) x (2 A +1): 

Sjv = I n) : —A < m, n < A 1 , (30) 



where m and n are integers. Here we make an assump- 
tion, well justified by numerical calculations, that the 
procedure produces convergent result as the number of 
plane waves in the expansion (|29|l is increased. Note that 
the choice of subset Eat of reciprocal lattice vectors Q is 
crucial, since the matrix elements decrease only as 1/|Q| 
for large momenta, and if the matrix in reciprocal space 
was truncated differently, the eigenstates could end up 
corresponding to a different self-adjoint extension. Af- 
ter substitution into Hamiltonian Tip x the coefficients 
= ( w Q' w q) can be found by diagonalization of the 
Fourier-transformed equations 

i 

#o(k + Q)tf£+ E F(Q - Q') = > (31) 

Q'eSjv 

where the summation is performed over all reciprocal vec- 
tors Q' e Sat except Q' = Q, the zero field Hamiltonian 
#o(p) is 

i?o(p) = cr 3 p x + aipy , 
and {mhvBz)~ 1 Q 2 V{Cl) is given by 

I fit ( e -iQ-B.B _ e -iQ RA) 2 Q^ e -iQ R B ) 



Using identities V(-Q) = V(Q) = -<r 2 V(Q)a 2 and 
a 2 H a 2 = —H Q it is easy to verify that if is the 
eigenfunction of l|31|) with momentum k and energy E 
then 

*Q = CT 2*Q> Q G Ejv 

is also a solution of (|31|l with momentum k and energy 
— E. In real space the new wavefunction corresponds to 
\J/'(r) = CP*(r) = a 2 ^(-r). Similarly, from identity 
V(Q) = — e 2< 3' R/1 V(— Q) we can obtain another eigen- 
state 

with with energy —E but inverted Bloch momentum 
— k. In real space the new wavefunction is described by 
= ^(2Ha — r). In fact using the properties of matrix 
elements V it is straightforward to show that all transfor- 
mations of Table II I II from Appendix that do not involve 
factors e~ lS ^ are exact symmetries of the Hamiltonian 
for any finite A. According to our assumption the 
procedure converges to a solution, which corresponds to 
some self-adjoint extension described by certain #1(2) • By 
inspection, we find from Table II I II from Appendix that 
the only choice of #1(2) ensuring the above symmetries 
are satisfied for arbitrary finite A is #1(2) = ±7r/4. 

We conclude by demonstrating that this self-adjoint 
extension has a single nodal point with zero energy at 
k = 0. Since Hamiltonian Hft does not possess a small 
parameter, and both the scalar and vector potentials v 
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FIG. 5: Left panel: the energy bands of linearized Hamiltonian for two vortices per unit cell calculated by plane wave expansion 
(solid lines) and by OPW method with boundary conditions 8a = —n/4 and 6b = +7r/4 (dots). Center panel: symmetry 
operations relating eigenenergies at different k. Right panel: the first Brillouin zone with equivalent k vectors shown by circles. 
For every state with energy E at momentum k there is a state with the same momentum and energy —E. 



and a are singular near vortices, it is desirable to base 
the discussion on non-perturbative argument. Such an 
argument was suggested by Vishwanathii^, who made 
use of two properties of the energy spectrum at k = 0. 
First, all states are doubly degenerate at k — 0, second, 
for each state (k, E) there is another state (k, — E). He 
noted that the two properties are shared by the zero- 
field problem, which has a pair of states at zero energy 
and symmetrically located pairs of states at E n and E_ n , 
where n is the band index. Thus, the number of states 
is 2(2n max + 1), where n max is a heuristic parameter 
denoting the number of bands with positive energy after 
fictitious ultraviolet cut-off is introduced in order to make 
the spectrum bounded. According to Vishwanath, since 
the total number of states is preserved after the potentials 
v and a are turned on, the pairs of degenerate states of 
TLft must also be centered around zero. Otherwise the 
total number states in the system would have been a 
2(27i maa; ). 

The argument in the form just presented might seem 
ill-defined and lead to inconsistent results. Indeed, if 
applied blindly, it would falsely suggest that the spec- 
trum should also be gapless at the corner of the Bril- 
louin zone. Certainly, the number of states at the center 
and the corner of BZ must be the same for any reason- 
able high-energy cut-off and therefore should be equal to 
2(2n max + 1). Just as for the states at the center of BZ, 
the states at the corner are doubly degenerate, and for 
each pair of states there is a symmetric pair with the 
opposite energy. Thus, one might conclude, there must 
be a pair of states precisely at zero energy not only for 
the Hft but also for the zero-field problem! Yet, the 
spectrum at the corner of BZ is gapped. 

Clearly, the difficulties one encounters by applying the 
above argument literally are related to the high-energy 
cut-off and have a technical character. We show below 



that Vishwanath's construction** can be made precise 
in the framework of the self-adjoint extensions adopted 
in this article. We have seen that the plane wave ex- 
pansion corresponds to the extension characterized by 
#i(2) = ±7r/4 if the set of reciprocal wave vectors of the 
basis forms a square grid (|30|l as N — > oo. At each finite 
N, there are (2N + l) 2 reciprocal lattice vectors Q con- 
tained in the grid £ jv and each Q is associated with two 
basis vectors (1,0) and (0,1). Thus the dimensionality 
of the basis is 2(2N + l) 2 . Since even for finite N the 
solutions at k = appear as doublets at energy E and 
another doublet at —E, the argument of Vishwanath can 
be directly applied, and one of the doublets should oc- 
cupy E = 0, otherwise the number of states would have 
been a multiple of 4 rather than 2(2N + l) 2 . 

It is interesting to observe why the same argument 
cannot be applied to the eigenfunctions of TLft at the 
corner of the Brillouin zone, and the spectrum remains 
gapped. Although at k = tt the spectrum is degenerate 
in the limit of infinite grid of reciprocal lattice vectors, for 
finite N the degeneracy is only approximate. As N is in- 
creased, the low energy bands become nearly degenerate, 
however at high enough energies the spectrum at k = tt 
remains corrupted. Transformation <?2 = D(CPI), which 
ensured the degeneracy of the spectrum at k = due to 
5252 = — 1 even for finite N does not have the same ef- 
fect on the states at k = tt. Certainly, after applying 
92 to an eigenstate xp„ with momentum tt and energy 
E, one still obtains an eigenstate ip-^ = c^Vv with the 
same energy, but with momentum — 7r. The situation 
is shown schematically in Fig. Although for infinite 
N the points tt and — tt are equivalent, for any finite N 
the wavefunctions describing the two states contain dif- 
ferent set of reciprocal lattice vectors, and therefore are 
distinct. In the analysis of the degeneracy of the states 
the two values of k should be considered as independent. 
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FIG. 6: Schematic representation of the spectrum at finite 
N. The small gaps between given two neighboring bands at 
the corner of Brillouin zone decrease to zero in the limit of 
infinite N. 



if Ri is an A vortex, and 
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e -r /a 
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if Hi belongs to class B. The wavefunctions that have the 
same asymptotics at the vortices and Bloch periodicity 
can easily then be constructed as 



XqO) 



R 



e^ X (r + R) 



The wavefunctions are orthogonal to the singular parts 
of asymptotics (llOfl . The appropriate basis for expan- 
sion of the eigenfunctions of the Hamiltonian TIft with 
boundary conditions {9i} is 



Thus the argument on zero-modes cannot be applied at 
the corner of Brillouin zone and for #1(2) = ±7r/4 there 
is a single Dirac node at P-point k = 0. 



V. NUMERICAL ANALYSIS 

As we have shown, the straightforward expansion of 
wavefunctions V'k(r) in the plane wave basis 



and 



iq-r 



where q = k + Q, Q e S is problematic since the ma- 
trix elements of the Hamiltonian decrease only as 1/Q 
because of the singular nature of the perturbation po- 
tentials at vortex locations. Consequently the eigenval- 
ues of the linearized Hamiltonian depend on the way the 
matrix (q'\H\q) is truncated and consistent ultraviolet 
cut-offs are parameterized by parameters 9 specifying the 
boundary conditions at each vortex. For example, the 
symmetric choice of the set Q €E Sat results in a choice 
&i(2) = ±7r/4. To enforce the desired boundary condi- 
tions we use a procedure similar to orthogonalized plane 
waves method used in the standard band theory. Instead 
of expanding the wavefunctions ip(r) in plane wave basis 
we define 



2 / 2 
-r I a 



sniffed 



ar 



e-*+ - 1 
ie'^ - 1 , 



cos tie 2 



cos tf> 



-i<f> + 1 



(32) 



E 



and 



where cutoff a is introduced for numerical convenience. 
The coefficients and /3 q are chosen so that the basis 
functions are orthogonal to x l ( r )- I n the limit a — > 
the procedure ensures that the modes with asymptotics 
orthogonal to those given by 9 are projected out. Be- 
cause of the exponential factors introduced in x( r ) the 
largest reciprocal vector Q should be of order 1/a or 
larger, so that there is appreciable region around a vor- 
tex 1/Qmax < r < a where the wavefunctions have the 
desired l/y/r asymptotics but still are not affected by 
the exponential factors. We typically used a = 0.051 and 
increased the number of basis elements until the con- 
vergence is achieved. Further decrease of the damping 
parameter a affects the results very weakly. Examples of 
the energy bands obtained by this method for two vor- 
tices per unit cell for different boundary conditions are 
shown in Figs. |S] and [7| The first figure shows com- 
parison of the energy spectra for the case 9\ = — 7r/4, 
9-2 = 7r/4 calculated by using the technique just described 
and straightforward plane wave expansion. The two are 
identical within numerical precision, confirming yet again 
that simple plane wave expansion generates a solution 
described by 0U2) = ±7r/4. 

Quite a different result is obtained for boundary con- 
ditions 9a = 0, 9b = 0. The spectrum in this case is 
gapped. The symmetry analysis can be applied to de- 
scribe the spectrum just as in the previous section. Rel- 
evant symmetry transformations and the star of equiva- 
lent symmetry points are shown in Fig. The features 
of the numerical results fully agree with the results of 
numerical calculation: for example the energy bands are 
symmetric under transformation {E, k) — > {— E, k), and 
the segments XM and YG are equivalent. 

As was pointed out, for general boundary condition 
even the symmetry (-E^k) (— E, k) of the single node 
energy spectrum is absent. 
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FIG. 7: Left panel: the energy bands of linearized Hamiltonian for two vortices per unit cell and boundary conditions 9a = 
Ob = 0. The spectrum obtained by OPW method is shown by solid lines, the dots represent the spectrum of Hamiltonian 1341 
with mass M(r) chosen as Gaussian Mo exp(— r 2 /£ 2 ), centered around each vortex with parameters Mo = 100//, £ = 0.05L 
Small discrepancies between the results obtained by two techniques for high energy bands are due to finite size of the core size 
£/l used in the numerical calculation. Center panel: symmetry operations relating eigenenergies at different k. Right panel: 
the first Brillouin zone with equivalent k vectors shown by circles. For every state with energy E at momentum k there is a 
state with the same momentum and energy —E. 



VI. RELATION OF THE SINGULAR 
LINEARIZED HAMILTONIAN TO THE 
NON-LINEARIZED, REGULAR PROBLEM 

As explained in previous sections, the spectrum of lin- 
earized Hamiltonian, in which vortices appear as point- 
like defects, depends on the boundary conditions imposed 
at vortex locations. These boundary condition are in 
turn determined by the self-adjoint extension of the lin- 
earized model. In principle, in order to establish which 
particular self-adjoint extension should be used, one faces 
solving the original, fully self-consistent BdG equations, 
a project beyond the scope of the present paper. 

Still, in the simplest case of a <i-wave superconductor 
with no additional forms of ordering in the interior of 
vortex cores, the physics of nodal quasiparticles should 
be adequately described by the Hamiltonian of Eq. @ 
and we should be able to say something definite about 
the meaning of different self-adjoint extensions. A useful 
perspective on various choices of boundary conditions, 
i.e. various choices of the parameter 9, can be gained by 
noticing that the main feature of the self-consistent solu- 
tion is to suppress the quasiparticle wavefunctions inside 
vortex cores. Furthermore, such fully self-consistent solu- 
tion does not require any additional boundary conditions 
at vortex locations, since non-linear terms, which grow 
strongly near vortices, and the self-consistency condition 
conspire to regularize the wavefunctions inside the core. 

How can such behavior be emulated on the level of 
a Dirac Hamiltonian without explicitly enforcing the 
boundary conditions? Unlike particles described by a 
Shcrodinger Hamiltonian, Dirac fermions cannot be pre- 
vented from penetrating vortex cores by a strong scalar 



potential barrier at vortex locations - if such a barrier 
is imposed it leads to the Klein parados^. The correct 
procedure which ensures suppression of the Dirac spinor 
amplitude in a particular region of space requires that 
the mass of the particle be treated as a function of po- 
sition M(r), so that in the prohibited region such mass 
becomes very large^. If the mass of the particle in the 
interior of forbidden region is set to Mo, then the spinor 
wavefunctions experience an exponential suppression in 
this same region, with a penetration length ~ 1/Mq. In 
order to regularize the FT Hamiltonian by requiring sup- 
pression of the nodal quasiparticle wavefunctions inside 
vortex cores, we are thus led to introduce a short-ranged 
mass-like potential <72-M(r) that vanishes at distances 
larger than the core size ~ £. Inside the core, the ab- 
solute value of mass Mq should be chosen to be much 
larger than l/£. 

Even in the zero-field problem such mass term breaks 
time-reversal symmetry: <72-M(r) changes sign under the 
time-reversal operation, and one might expect that in 
such terms in general lead to opening of a gap in the 
spectrum. The detailed analysis shows that FT equations 
augmented by mass potentials 

—TL' ft = (Px + a x )<T 3 + a A a 1 (p y + a y ) +v x + a 2 M(r) 

Vp 

(34) 

are generally not invariant under symmetry operations 
listed in Table II I II from Appendix. The transformation 
properties of the mass term itself M(r) are shown in Ta- 
ble [H] Clearly, no choice of nonzero mass term M(r) 
is invariant under all symmetry operations, just as no 
choice of O^s preserves all the symmetries of Hpp. For 
M(r) possessing certain specific symmetries, some of the 
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1 I P PI 

M(r) -M(Zr) -M(Pr) M(sr) 
C I CI I CP I CPI 
-M(r) M(7r) M(Pr) -M(sr) 

TABLE II: Transformation of the position-dependent mass 
term M(r) under operations commuting with Hft- Only 
for M(r) of special symmetries the new Hamiltonian H' FT = 
Hft + 02VfM(y) remains invariant upon the action of these 
operations. 

transformations do commute with the Hamiltonian. For 
example, if the mass does not change under inversion 
about a vortex and has the same sign for all vortices: 

M(Jr) = M(r) = M(sr), 

then the transformations commuting with the Hamilto- 
nian coincide with those listed in Fig.[7J Thus the Hamil- 
tonian with such mass terms is equivalent to the original 
FT Hamiltonian with 8\ = 82 = and with mass equal to 
zero. The spectrum of Ti' FT calculated by the standard 
plane wave expansion is shown in Fig. and agrees with 
the numerical results based on the OPW method. Simi- 
larly, a straightforward calculation shows that the choice 
of 

M(r) = M(Jr) = -M(Pr) 

corresponds to 81 = Q; 62 = n/2, while 

M(r) = -M(Jr) = -M(Pr) 

describes 6\ = —82 = tt/4. The latter choice requires 
that M(r) has a p-wave symmetry and vanishes at least 
along two directions around each vortex. This choice 
also illustrates the following important point: while the 
mass term always breaks time-reversal symmetry, in the 
8\ = —82 = 7r/4 case it does so only locally while leav- 
ing this symmetry intact on average since J d 2 rM(r) = 0. 
Here the integration is over the small region within which 
M(r) differs appreciably from zero around a single vor- 
tex. 

We will conclude this section with two observations. 
First, if the overall amplitude of the mass term is in- 
creased continuously, then the gap in the spectrum will 
also increase, and it is not a priori clear that in the 
limit of extremely small magnetic field the magnitude 
of the gap induced by the mass term is proportional to 
l/l rather than l/l 2 . To answer that question we note 
that if the absolute value of the mass M (r) inside vortex 
core is denoted as Mq, then the quasiparticles penetrate 
into the vortex core to distances of order 1/Mq. Since 
the distance should be much smaller than the size of the 
vortex core £, we find that Mq 3> Z/£. The scaling limit, 
therefore is obtained as a limit of infinite Mq. 

Our analysis does not allow us to uniquely determine 
which self-adjoint extension (specified by 8\^ or by sym- 
metry of the mass terms) is realized in cuprates, as the 



answer depends significantly on short range physics and 
can be determined only through the detailed analysis of 
the spectrum and the wavefunctions of self-consistent so- 
lution. However, given the extension of the Hamilto- 
nian describing quasiparticles near node l, the symmetry 
properties of the full BdG Hamiltonian determine unam- 
biguously the self-adjoint extensions of linearized Hamil- 
tonian near nodes 1, 2 and 2. 

Symmetry operations of the linearized Hamiltonian 
listed in ( Table ITTT1 in the Appendix) also commute with 
the full BdG Hamiltonian Ti defined by JSJ, where now 
these symmetry operations are applied to the full wave- 
function 'J. One should note, however, that these op- 
eration have entirely different meaning as Ti describes 
simultaneously all four nodes and most of the operations 
relate the eigenstates belonging to different nodes. For 
example, the action of the inversion operator on an eigen- 
function \&(r) of Ti will now be understood as , J , (— r): 

¥(r) « e lkpX ^{v) -» e- ikFX iP{-r) ; 

the reader should compare this to the "inversion" oper- 
ation of Eq. H25fl . Thus, it transforms the state near 
node 1 into a state with Fourier components localized in 
the vicinity of 1. Similarly, symmetry operations involv- 
ing reflections m x and m y applied to quasiparticle wave- 
functions with Fourier components localized near node 1, 
generate new wavefunctions near nodes 2 and 2. 

For simplicity, we will focus on transformations gener- 
ated by P, I, and C, which relate nodes 1 and 1 only. 
Acting on a wavefunction ijj(r) with crystal momentum 
k and energy E, they generate seven new states shown 
schematically in Fig. [H] 
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FIG. 8: Symmetry of the spectrum for the full BdG Hamil- 
tonian Ti. 

Note that transformations PI, CP, and CI relate 
wavefunctions near the same node. They ensure that 
if ^(r) is an eigenstate of Hamiltonian J2J) with energy E 
and momentum k then the spectrum of Ti also has states 
characterized by (k, —E), (k + tt,E), and (k + 7r, — E). 

One can easily verify that there are only two self- 
adjoint extensions compatible with PI, CI and CP: 
01 = &2 = and 8 \ = 82 = n/2. The two choices are 
essentially equivalent being related to each other by re- 
flection m x . The resulting single-node energy spectrum 
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is shown in Q. It is characterized by strongly disper- 
sive bands at low energies, except for the lowest gapped 
narrow bands at energies w ±2Hvf/1- 

All other boundary conditions break one of the three 
symmetries. In particular, if the spectrum is gapless at 
kf then it should also be gapless at momentum kf + tv. 
This is clearly not the case for the choice 0\ = 7r/4, 
62 = — 7r/4 shown in Fig. |SJ To restore the full sym- 
metry of the non- linearized BdG Hamiltonian J3J, one 
has to demand that the linearized version is a union of 
two independent self-adjoint extensions: Q\ = —62 = 7r/4 
and 9\ = —62 = — 7r/4. Such union is a combination 
of the solution found in Ref. and its companion so- 
lution with the A and B labels interchanged (DOS, of 
course, remains unaffected apart from a mere factor of 
two). Consequently, the single- node linearization cor- 
responds to two independent problems described by the 
same Hamiltonian operator, but with different bound- 
ary conditions near vortices. Under PI operation, the 
Hilbert spaces of the two Hamiltonians are interchanged 
and the full symmetry is thereby restored. Observe that 
this is a deeply non-perturbative result: such a union de- 
scribes eight two-component Dirac fermions and therefore 
16 zero energy states at BdG nodes, in contrast to only 
four massless Dirac fermions of the H = case. Recent 
exact symmetry arguments on a tight-binding lattice 7 - , 
based on a form of the index theorem, strongly support 
the identification of the above union as the proper self- 
adjoint extension of the continuum linearized problem 
with gapless spectrum. 



VII. LINEARIZATION OF d x 2_ y 2 
HAMILTONIAN 

In the companion paper—, we consider properties of a 
tight-binding lattice d-wave superconductor in the mixed 
state. The simplest representative of such systems is 
characterized by the d x 2_ y 2 symmetry of the order pa- 
rameter rather than the d xy case considered in previ- 
ous sections. In order to directly compare the linearized 
model and the tight-binding lattice model, we use this 
section to consider the linearized Hamiltonian with the 
"clover" of the d-wave gap function rotated by 45 degrees 
with respect to the vortex lattice. After following the 
standard stepaiS<2I, we find that the linearized Hamilto- 
nian in this case is given by 



rrlin 

n x 2_ y 2 



n r + n„ 



n„ - dl 



v F - 



V2 



-OS + VA- 



V2 



-<7i + Vf- 



V2 



(35) 

where = pi + a$ is the generalized momentum. In 
the rest of the section, we consider the properties of the 
above Hamiltonian for precisely the same orientation and 
position of the vortex lattice as before, with the unit 
cell shown in Fig. ^ Therefore, the potentials v and a 
remain exactly the same, and the only difference from the 
original Hamiltonian TLft is the direction of the d-wave 
nodes relative to the vortex lattice. 



As before, we will confine ourselves to the isotropic case 
vp = i'a- Note that in this case, if one neglects the super- 
fluid velocity terms and retains only the vector potential 
a, Hamiltonians 2 and Hft = H l xy are related by a 

unitary SU(2) transformation U = exp(icr y 7r/8). There- 
fore, both are expected to have the same density of states. 
Moreover, since this unitary transformation is global, the 
dispersions E n (k) of the two cases will also be identical. 

The above invariance is a consequence of the combina- 
tion of our Dirac particles effectively being at k = and 
the axially-symmetric character of the free Dirac prob- 
lem. Since there is no preferred direction for the free 
Dirac problem, different orientation of the vortex lattice 
cannot change the spectrum. This invariance, of course, 
holds for any linearized Hamiltonian with an arbitrary 
direction of the d-wave "clover" . The unitary transfor- 
mation in this general case, relating the Hamiltonian to 
Hft is U = exp(io~ y a/2), where a is the direction of the 
node relative to x-axis. 

Of course, this invariance is violated by the superfluid 
velocity terms and finite anisotropy. In the latter case, 
there is a preferred direction since the free Dirac dis- 
persion defines an elongated cone. Similarly, the terms 
containing the superfluid velocity v carry the informa- 
tion about a particular direction of the node, where the 
linearization is performed, and therefore break the rota- 
tional symmetry, even when a is set to zero. Numeri- 
cally, the effects of v turn out to be rather small at the 
lowest band energies, and become progressively more im- 
portant for higher bands. Fig. El shows the spectrum 



of cL 



superconductor obtained numerically using a 



simple expansion in plane waves. For the lowest bands, 
it is very similar to the spectrum of dry-superconductor 
discussed in detail in the earlier sections (see Fig. 0). In 
both d xy are d x 2_ y 2 cases, for the lowest band the result 
is very close to that of the reduced Hamiltonian, which 
contains only the vector potential a - the effect of the 
scalar potential becomes pronounced only at higher en- 
ergies. The spectrum of the Hamiltonian with only the 
scalar potential v included and the vector potential set 
to zero by hand, is quite different - at low energies it is 
a tiny perturbation of the free Dirac dispersion. This is 
in stark contrast with the result of the so-called "Volovik 
approximation . 

The allowed class of the boundary conditions near vor- 
tices for H l ^2_ y 2 can be obtained directly from the bound- 
ary conditions of the d xy case (1221 and (123(1 by noting that 
in the coordinate system (x',y') rotated by 45 degrees 
with respect to (x,y), H^P 2 takes on the form identi- 
cal to Hft , except for the vortex lattice being rotated by 
7r/4 relative to what it used to be in the d xy case. Still, 
when deriving the asymptotics of the wavefunctions near 
a vortex, only the singular terms due to that single vor- 
tex contribute - different geometry of the lattice therefore 
does not directly affect the boundary conditions. All we 
need to do is to return to the original coordinate system 
by performing the rotation <j> — > — 7r/4in (|22[) and i|23ll . 

The symmetry properties of the Hamiltonian H l *™_ 2 
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FIG. 9: The quasiparticle spectrum of the mixed state ob- 
tained by an expansion of the wavefunctions in the plane wave 
basis. The dispersion of the zero magnetic field problem is 
shown in red squares, the dispersion of a d x 2_ y 2 superconduc- 
tor in the presence of a vortex lattice is shown by a solid black 
line. The green diamonds (blue triangles) correspond to the 
artificial Hamiltonian, which is obtained by setting the scalar 
potential v (vector potential a) in Hft to zero. Note that the 
role of the scalar potential at low energies is small, and that 
makes the spectrum of ffj,™ i rather similar to HJy = Hft 
shown in Fig. 

can be derived rather straightforwardly, since v and a 
are exactly the same as before, and we can use their 
properties from Table ^ The symmetry properties of 
the Hamiltonian, which do not involve mirror symme- 
tries, remain exactly the same. Namely, transformations 
J, P, C and their combinations can be literally taken from 
Table IIIII The transformations involving mirror symme- 
tries arc modified, however. Moreover, in addition to the 
mirror planes along the x and y directions, now a new 
mirror symmetry plane appears along the diagonal of the 
unit cell. The appearance of these new elements of sym- 
metry is important for a precise form of the spectrum. 
It is easy to see from Figs. El and El that, although the 
spectra of the d xy and d x i_ y i are quite close, the latter is 
more symmetric due to the additional symmetry x <-> y. 

The main conclusions, nevertheless, remain the same. 
For example, the symmetry of the Hamiltonian operator 
requires that the spectrum is symmetric under transla- 
tions in momentum space by (tt/1,tt/1). In other words, 
if there is a node at k = 0, it should have been replicated 
at the corner of the unit cell. Clearly, the dispersion of 
Fig. El has only one gapless point at k = 0. The reason, 
just as before, is the necessity of imposing the boundary 
conditions which again violate the symmetries of what is 
formally the Hamiltonian operator. 



VIII. CONCLUSIONS 

The main results of this paper can be summarized 
as follows: the linearized BdG Hamiltonian describing 



nodal fermions in the mixed state of a d x 2_ y 2 supercon- 
ductor has to be complemented by a set of boundary 
conditions Ij22l23(l specifying the behavior of quasiparti- 
cle wavefunctions near vortices in order for the problem 
to be mathematically fully defined. The boundary con- 
ditions contain a single parameter 9 € [0, n) for each 
vortex which cannot be found from within the linearized 
theory itself, but should be determined from the full non- 
linearized calculation. All the physics beyond lineariza- 
tion, such as the intervortex scattering and interference 
effects, particle-hole asymmetry, high energy and curva- 
ture terms, etc., is implicitly reflected in the effective 
linearized FT Hamiltonian only through these boundary 
conditions. Consequently, the linearized FT Hamiltonian 
actually describes a family of distinct self-adjoint exten- 
sions reflecting a variety of high-energy processes that 
might be taking place inside the vortex cores (magnetism, 
Mott insulator, charge density-wave, etc.). Such multi- 
tude of all possible short-range physics behaviors can be 
classified according to the set of 9 parameters, which gov- 
ern different "universality" classes of the nodal fermion 
quantum "criticality." The conventional Simon-Lee scal- 
ing function must be generalized to include the explicit 
dependence on 8's, as explained in the text. Once the 
boundary conditions are properly imposed according to 
(I22I23[1 . the solutions of the FT equations are fully de- 
termined and are independent of the assignment of A or 
B vortices, thereby restoring their invariance under arbi- 
trary singular gauge transformations. 

In earlier work 8 the boundary conditions were not ex- 
plicitly enforced but instead were "spontaneously" se- 
lected by the procedure used to diagonalize the Hamil- 
tonian. Such an expansion in the plane wave basis re- 
sults in 9u2) — ± 7r /4(T 7r /4), which was indeed found 
here to be the self-adjoint extension appropriate to the 
nodal gapless behavior of a d-wave superconductor in the 
limit of low magnetic fields. In general, however, the 
boundary conditions #1(2) should be found from the full 
self-consistent BdG Hamiltonian and serve as an external 
input to the linearized theory describing the bulk quasi- 
particle states. 
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and Profs. M. R. Zirnbauer and A. Altland for helpful 
correspondence and for sharing with us their unpublished 
results. This work was supported in part by the NSF 
grants DMR-0094981 and DMR-0531159. 



APPENDIX A: SUPERFLUID VELOCITY AND 
THE PHASE OF SUPERCONDUCTING ORDER 
PARAMETER 

1. Simple properties 

For the reference, we summarize the properties of su- 
pcrfluid velocity and the phase of the order parameter 
used in the text. Consider a periodic lattice of vortices 
with a basis, with N vortices located at within a 
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square unit cell of size I x I as shown in Fig. ^] The 




FIG. 10: Left panel: An example of a vortex lattice for 
the case of N=4 vortices per unit cell. Right panel: Two vor- 
tices per unit cell: integral l!A17fl calculated along the contour 
shown as a dashed line equals 2n. 



superfluid velocity is defined as 
tov = —fiVcf) — 



(Al) 



where </>(r) is the phase of the order parameter. The 
latter can be determined from the requirement that it 
acquires 2tt as r encircles a vortex: 



V x V0 = 27rz^(5(r-R) 



(A2) 



R 



where R denotes the position of the vortices. Since the 
superfluid velocity typically enters in combination mv we 
will often set m = 1, for compactness. 

Using the Fourier transforms, the superfluid velocity 
can be written^ as 



v(r) = iirh 



ark kxi 
(2tt) 2 A- 2 + k 2 



R 



-ik-R 



(A3) 



where A is the London penetration length. Let us write 
the position of a vortex as R = R ; + r, where index 
i = 1. . .N labels the vortices inside an arbitrarily cho- 
sen reference unit cell, and r = l(N x x + N y y) with in- 
teger N x , N y denotes lattice translation vectors. Then, 
after summation over r, the superfluid velocity can be 
expressed as 



v W = -p E 

Q 



{Qyi Qx) 

Q 2 + X- 



N 



iQRi 



(A4) 



N y are the 



where Q = 2ir(xN x +yN y )/l with integer N L 
reciprocal lattice vectors. Note that when X^> I, one can 
neglect A~ 2 in the denominator in all terms of the sum 
except for Q = 0, and the simplified result reads: 



v ( r ) = -p E 

Q/0 



(Qy, Qx) iQ, 

r>2 e 



N 

E< 

i=l 



-iQRi 



(A5) 



This simplification corresponds to replacing the magnetic 
field H(r) by its spatial average, as can be verified di- 
rectly by comparing V x v from (|A4I) and 1|A1(1 . In the 
rest of the appendix and in the main text, the limit A I 
is always implied, and A(r) denotes the vector poten- 
tial corresponding to the uniform applied magnetic field, 



without including small corrections due to the vortex lat- 
tice, which are smaller by a factor of l 2 /X 2 . 

The superfluid velocity v can be expressed in a closed 
form through Weierstrass elliptic zeta function (see also 
Ref. l55h . To derive this expression, let us start from 
(|A2(I . and write the phase <fi as 



(r) = Mr) 



y arctan 



R 



X 



X 



where <fio(r) is a continuous function to be determined 
later and the second term is a sum of the polar angles 
describing r with respect to the vortex location R. After 
denoting z = x + iy and substituting <fi into (|A1|) we 
obtain 



N 



Vy + iv x = 



1 



z - Zt 



(Ay + iA x ~) 



+ 2^1/^0 + ^00) , (A6) 

where the lattice vectors are r = t x + ir y and the loca- 
tion of the vortices within the reference unit cell are now 
encoded by Z% = X t + iYi. The sum over the lattice vec- 
tors t in the right hand side differs from the definition of 
Weierstrass zeta-function 



C(*) = - 

z 



E 



1 



(A7) 



only by a presence of the constant and linear in z 
terms, which are required to ensure the absolute con- 
vergence of the sum <|A7ll . These terms can be ab- 
sorbed into the "smooth" part of the phase by defining 
= 0o + Cqz + \C\z 2 , where are appropriately 
chosen constants. Thus, the superfluid velocity satisfies 
the following equation: 

%+-*=!e(c(*-^)-^ 1 v »* ,+<v "# 



(A8) 



Here, for convenience we chose the symmetric gauge: 



£ -A 



7TN 

2p (-¥.*) 



satisfies both 



Note now that w{x, y) — Vj,0 o + i\7 x 
Cauchy-Riemann conditions: d x (Rcw) = d y (lmw) is 
satisfied automatically, whereas the second condition 
d y (Rew) = —d x (Imw) is satisfied due to the requirement 
div v = and i|A8(l . Since 4>o (and consequently <j>' ) was 
chosen as a smooth part of the phase, it is therefore a 
finite analytic function of z in the entire complex plane, 
and by Liouville theorem must be equal to a constant. 
The remaining constant is fixed by the requirement for 
the spatial average of v to vanish (cf . IA5|) , and the final 
result, which is equivalent to l|A5|l . reads 



* N 
2 ^ 

i=l 



C(z - Zi) - irh- 



Zi 



(A9) 
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Note that the superfluid velocity v(r) is periodic in space Eq-n l|A12|) is the explicit expression for the phase <p(r) 
with the same unit cell as the unit cell of the vortex used in the main text and in the rest of the appendix, 
lattice. 



2. The phase of the order parameter 

Unlike the superfluid velocity components, neither the 
phase of the order parameter <fi(r), nor exp(i</>), and not 
even V</> are periodic, regardless of the gauge used. Since 
the periodicity of exp(i0) would have required that 
is itself periodic, it is sufficient to prove the statement 
for \7(f>: consider a contour integral <f V</> • dl along the 
boundaries of the unit cell. Since there are N vortices 
inside the contour, the integral must equal 2?tN. On 
the other hand, the assumption of V0 being periodic 
would have resulted in vanishing of the integral due to the 
cancellation between the contributions of the opposite 
edges. 

Although the phase cannot be made periodic, it is 
quasi-periodic and can be expressed in a closed form 
through Weierstrass sigma-function a(z). Since (|A1|I 
yields 



d v <f) + id x 4> = ^U(z-Z i ) + 



the difference cf>(r) — 4>(yq), where ro is an arbitrary ref- 
erence point, can be written as 



V</> • dl = Im / (d y (f> + id x (j))dz 

J ZQ 

' Z TtZAz — Zn) i 

tic + — I , (A10) 



= Im E ( f Z "cm<* 

.■ V Jz -Z, 



I 2 



where the contours in all integrals are assumed to be the 
same. For different contours, the equality holds modulo 
2ir. Since the Weierstrass sigma-function a(z) is defined 
according to a'(z)/a(z) = ((z), which implies 

| n l^L = j C,{z)dz (mod 2tt) 



O"0o) " Jz 

the phase <fi(r) can be written as 
0(r)-^(r o )^Arg g(z -y 



-Im I - zo 



J2 Z ^j ( m ° d 2?r ) • ( AU ) 



Since the phase <^(r) is defined up to an arbitrary con- 
stant, even for a fixed gauge of the vector potential A, it 
is convenient to choose the constant </>(ro) in such a way 
that 

4>(r) = Y^Atg[a(z-Zi)]+^lm i^zZ^\ ( mod 2?r ) ■ 

(A12) 



3. Two vortices per unit cell 

In general, the procedure of separation of the phase 
</>(r) into (f>A and 4>b is straightforward by using l|A12|l . 
In this section we illustrate it for the simplest case of 
N = 2 with two vortices located at Ha and Rb inside 
an arbitrary reference unit cell of size I x / (see Fig^l. 
The vector potential in the symmetric gauge is given by 

e a n i ^ 

-A = -p>\—y> x ) , 

c l l 

and the velocities va,b can be defined similarly to v: 



v A(B)(r) 



Zi^h X - (Qy, Qi) ,Qr -iQR A(; 
72 O 2 



(A13) 

Clearly, v A (b) automatically satisfy va + vb = 2v (see 
(|A3(l h Moreover, since Va(b)/% is formally given by the 
same Fourier expansion as v for N — 1 (cf. HA3|l ). we 
have 



V A(B) _ 1 



V0A(S) 



(A14) 



2 2 ™ w 2c 
where 4>a(b) satisfies the analogue of l|A2(l : 

Vx V0a(b) =27rz^<J(r-R A(B) -T) , (A15) 

T 

The closed form expressions for \~a(b) can be read off 
directly from the formulae of the previous section. After 
defining 

F(r) = Arg[a(x + iy)] 



Z(r)= (ImC(z) + 1 i,ReC(z)- 1 3- 



where cr(z) and £(z) are Weierstrass functions with pe- 
riods the superfluid velocities are simply va(b) — 
Z(r — Ha(b))- The phases 4>a(b) are given by 

^ (s) (r)=F(r-R A(s) ) + ^A(R A(B) )-r ; (A16) 

by construction they obey the condition 4>a (r) + 4>b (r) = 
<Kr). 

In our discussion of the quasiparticle spectrum, we use 
the symmetry properties of the phase difference 5<j) = 
cf>A(r) — 0s(r). Let us show that after a translation by the 
primitive lattice vector Ix, the phase difference acquires 
a 7r-shift. Consider the following integral along a straight 
line connecting r and r + Ix: 

J(r) = / {v A -v B )-d\ . 
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Functions ~va( B ) are periodic in unit cell and therefore 
I(x,y) does not depend on x: 

7(r)= ! [v A {^y)-v B {Um ■ 
Jo 

Next, we note that the contour integral along the path 
shown in Fig^Jis a multiple of 2ir: 

<j> (va - v s ) • d\ = j> (V0a - V0 B ) ■ d\ = 2ir(n A - n B ) , 

(A17) 

where the n A — n B is the difference between the number 
of A and B vortices inside the contour. Since the integral 
along the vertical edges vanishes due to the periodicity 
of superfluid velocities va( B ), we find that 

I( x >Vl) - I ( X )V2) = 2vrn . 

Choosing y\ and y2 so that the horizontal segments 
are symmetrically located around vortex A as shown 
in Fie llOl we have I(x,y%) = —I{x,y2) which yields 
I(x,yi) = 7r, and consequently, 



(r + Ix) — 5<fi(r) = 7r (mod 2ir) 



(A18) 



Combining this result with a similar identity for the 
translations in y-direction: 

S(j)(r + ly) — 8<j){r) = n (mod 2n) . 

we find that exp(i<5</>(r)) can be written as a product of 
exp(i7r(a; + y)/l) and a periodic function with a unit cell 
I X /. 

Other useful identities involving 5<j> — 4>a — 4> B are 



S4>(x,y)+8<f>(l/2-x,y)=ir/2 , 
5cj)(x, y) + 5cj)(x,l/2 - y) = -vr/2 

W( x , y) + 5<t>(-x, -y) = . 



(A19) 
(A20) 
(A21) 




= x 



FIG. 11: Elliptic coordinate system. Each line of constant 
2 is an ellipse. At z = the ellipse is reduced to a segment 
(—a, a) of the real axis. 



the technique to the case of two-component Dirac parti- 
cles. 

In what follows, we will extensively use Mathieu func- 
tions. Therefore, we start by reminding the reader of 
some of their essential properties: the details can be 
found in Refs. Fll[72l andFil Mathieu functions appear 
naturally as solutions of problems possessing elliptic sym- 
metry and expressed in elliptical coordinates. The most 
common examples are Laplace and Helmholtz equations 
with boundary conditions specified on an ellipse. The 
elliptic coordinates are defined according to 



x = a cosh z cos < 
y = a sinh z sin <f. 



(Bl) 
(B2) 



where a is a constant. Points with z = lie on a segment 
(—a, a) traversed twice when <fi is increased from to 27n 
first from +a to —a, and then backwards. For finite z 
the contours z — const are ellipses (see Figlll|): 
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APPENDIX B: TWO AHARONOV-BOHM 
FLUXES - EXACT SOLUTIONS 

1. Elliptic coordinates and Mathieu functions 

In this appendix, we analyze the properties of a two- 
component Dirac particle moving in the field of two 
Aharonov-Bohm fluxes. We will prove that the exact 
solution of the problem is not unique and depends on 
boundary conditions imposed near the fluxes. To achieve 
this goal, we first derive a method allowing us to con- 
struct exact solutions for the problem of a Schrddinger 
particle in the field of two Aharonov-Bohm solenoids, 
which carry arbitrary fractional fluxes a^o and f3(f>o , with 
4>o = hc/e. We illustrate our procedure by analyz- 
ing several simple, physically interesting examples with 
a, (3 — ±i. After solving several model problems for a 
Schrddinger particle in the field of two fluxes, we apply 



= 1 



a 2 cosh z a 2 sinh z 

and angular coordinate 4> is used to unambiguously spec- 
ify the position of the point on the ellipse of constant z. 
Each ellipse is described by a major and minor semi-axes 
of lengths a cosh z and a sinh z respectively. The linear 
eccentricity defined as the distance from the center to 
either focus, is the same for the entire family of ellipses 
and equals a. At large distances r a, the contours 
of constant z reduce to circles, and the elliptic coordi- 
nates effectively reduce to the ordinary polar coordinates 
(lnp, 0). Both the Poisson and Helmholtz equations are 
separable in these coordinates - this is why the elliptic 
coordinate system is useful in various problems of math- 
ematical physics in the first place. The Laplacian in el- 
liptic coordinates assumes the following form: 



dx 2 



dy 2 



d 2 



a 2 (cosh 2z — cos 2(f>) \dz 



d 2 
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9' = D(g)9 
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k 


01(2) 


m x 
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-E, 


77 — k 
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Pm x 
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-E, 




-02(1) 


Pm y 


ffi*(5r) 


E, 


( — k X , fcy) 


7T/2 + 02(1) 


PI 


e- <5 ^( fl r) 


E, 


77 + k 


02(1) 



(a) 





D(Cg)^ = *' 


E' k' 


01(2) 


c 


e - i ^ (T2 "$( r ) 


E, 77 - k 


7r/2 + 0i (2) 


Cm x 


a 3 *(gr) 


73, (k X ,—ky) 


01(2) 


Cm y 


a 3 *(<?r) 


— E, (—k x ,k y ) 


TT/2 - 1(2) 


CI 




-E, 77 + k 


-01(2) 


CP 


a 2 *(<?r) 


-75, k 


-02(1) 


CPm x 




-E, 77 + (— fcs, k y ) 


TT/2 - 2( i) 


CPm y 




E, 77 + (fez, -fcj,) 


02(1) 


CPI 


a 2 *(<?r) 


73, k 


7T/2 + 02(l) 
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TABLE III: Transformation of the wavefunctions and the self-adjoint extensions under various symmetry operations. The right 
column contains operations involving complex conjugation. 



and the Hclmholtz equation V 2 / + k 2 f ~ becomes 



d 2 d 2 



dz 



2 1.2 



a 2 k 



+ -^72 f-\ (cosh 2z- cos 2</>)/ = 



After separation of variables via / = F(z)G(<fi), one ob- 
tains 



1 d 2 F 1 d 2 a 2 k 



2 7,2 



F{z) dz 2 G(0) d<t> 2 2 
which results in 

1 d 2 F a 2 k 2 



■(cosh2z — cos 20) = 



F(z) dz 2 2 
1 d 2 a 2 k 2 



G{4>) d4> 2 



cosh2z = ^ , (B3) 
cos 20= -A , (B4) 



where A is the separation constant. The conventional 
form of these equations, which of course must be solved 
simultaneously for the same constant A, is 

d 2 G 

— + (A-2q cos 20)G(0) = , (B5) 



d 2 F 
dz 2 



-(A- 2qcosh2z)F(z) = , (B6) 



where q is a parameter defined as q — a k : /4. 



a. Mathieu functions 

Since the solutions must be single-valued functions, 
G(0) must be 27r-periodic. This condition restricts A to 
a discrete set of so-called characteristic values, which is 
traditionally written as a union of two subsets represent- 
ing solutions that are even or odd under transformation 



ao, cii, 0,2, &3, ... (even solutions) 
bi, fe, ^3, ■■■ (odd solutions) 



While working with the elliptic coordinates, it is useful to 
keep in mind similarities and distinctions from the polar 
coordinate system. In the latter, the equation for the 
angular component assumes a form G" + AG = 0, and 
the requirement of 27r-periodicity restricts A to a set of 
m 2 with integer m, irrespectively of k. In the elliptic 
case, this no longer holds, and the characteristic values 
a.j, bj depend on q = a 2 k 2 /A. 

For A — a n the periodic solution of the equation 
for G(0) is an even function of 0, which is denoted as 
ce„(0, q), n — 0,1,2,.... The odd solutions, occurring 
for A = 6„, are denoted as se„(0, q). Functions ce n and 
se n are known as Mathieu functions, the notation is to 
remind us that they reduce to cos(n0) and sin(n0) re- 
spectively in the limit q — > 0. 

The equation for G is of second order, so for a given 
A = a n (or A = b n ) there is a second independent so- 
lution, in addition to ce n (4>,q) (or se„(0, q)). However, 
the second solution of the angular equation is never pe- 
riodic, and therefore it is of limited, if any, importance 
in applications. The characteristic values a n , b n , as well 
as the Mathieu functions ce n and se n are tabulated and 
well studied, see for example Ref. |^. Finally, we note 
that Mathieu functions {ce„(0, q), se„(0, q)} form a com- 
plete set in the interval (0, 27r) and satisfy the following 
orthonormality relations: 



2tt 



ce n (0, <?)ce m (0, q) d(f> = S, n 



se„(0, g)se m (0, q) dej) = 5 n 



2tv 



ce n ((j), q)se m {(t>, q)d(/) = 



b. Modified Mathieu functions 



(B7) 



(B8) 



(B9) 



Now we turn to the equation (|B6|) describing the "ra- 
dial" part of the solutions. Since (|B6|I and (|B5|I must 
be solved simultaneously for the same value of param- 
eter A, we are interested in solutions of i|B6(l only for 
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A 


1st solution 


f(0,q) 


9»/(0,«) 


a n 

bn 


Ce„(z,q) oc Afc (1) (:z,g) oc Je n (z,q) 
Se„(z,q) oc Ms w (z,q) oc Jo n (q,z) 








A 


2nd solution 


/(0,9) 


9*/(0,«) 


a n 
b n 


Fey n (z,q) oc Mc (2) (z,q) oc Ne n (z,q) 
Gey n (z,q) oc A'/s (2) (z,g) oc No n (z,q) 


7^0 


7^0 



TABLE IV: Modified Mathieu functions: notation and values 
at z = 0. 



Je 2n (z,q) s 


- P2n7 e ~ 


z/2 


sin(yge z + f) 


Ne 2 n(z,q) s 


- —pinker 


z/2 


cos(^e z + f ) 


Je 2n+ i(2,g) R 


" -P2n + l7 e ~ 


z/2 


cos( 1 /ge z + f ) 


iVe2„+i(z,g) s 


d — P2n + l7 e ~ 


z/2 


sin(yqe z + f ) 


Jo 2n (z,q) s 


- S 2 n7 e ~ 


z/2 


sin(^e z + f ) 


No2 n (z,q) s 


- — S2n7 e ~ 


z/2 


cos(y?e z + |) 


Jo 2 n+l(«,3) s 


d — S2„ + l7e" 


z/2 


cos(^e z + f ) 


ATo 2 ,i+l (z, ?) i= 


- — «2n + l7 e ~ 


z/2 


sin(yqe z + f ) 



TABLE V: Asymptotic expressions for Modified Mathieu 
functions for 2 > 1 in terms of 7 = 2 1/ ' 2 (7r 2 q)~ 1 ' /4 and the 
standard coefficients p n and s n (see Ref. r7Jl . 

A = a n or b n . These solutions are known as Modified 
Mathieu functions. For every value of parameter A there 
are two linearly independent solutions: for A = a n these 
solutions are denoted as Je n (z, q) and Ne n {z, q), and for 
A = b n the solutions are Jo n (z, q) and No n (z, q). These 
functions take place of Bessel J m (r) and Neumann N m (r) 
functions in the more familiar case of polar coordinates; 
the letters "e" and "o" stand for "odd" and "even" re- 
spectively 

Unfortunately the notation used for Modified Math- 
ieu functions is not standardized, and for convenience we 
provide the Table HVI which relates notations used by dif- 
ferent authors. The functions of the first kind, Je n (z,q), 
Jo(z,q), are proportional to the Mathieu functions of 
imaginary argument ce n (iz,q) and se(iz, q). In our sub- 
sequent analysis, we will also use the properties of the 
Modified Mathieu functions at z = and z ^> 1, which 
are summarized in the right two columns of Table ITVl and 
Table respectively. In the most general case, the so- 
lution of the Hclmholtz equation can be written in the 
following form: 

/ = ^ ce n (<f>, q)(a n J e n {z, q) + (3 n Ne n (z,q)j 

n=0 

+ se w (0, q) (j n Jo n (z, q) + S n No n (z, q) j , (B10) 

71=1 

where constants a n , n ,j n ,S n depend on the boundary 
conditions. In addition to the external boundary condi- 
tions, additional attention should be paid to the segment 
(—a, a) corresponding to z = 0: ordinarily the solutions 
and their derivatives are continuous across the segment 
providing a restriction on the coefficients cv n . . .5 n . We, 



however, will be also interested in solutions that experi- 
ence discontinuity across the line, which will provide a 
different set of constraints on the coefficients in (IB10|) . 

c. Model Schrodinger problem: a particle inside an elliptic 

box 

We begin with the wave equation describing a parti- 
cle inside an impenetrable box of elliptical shape. The 
wavefunction ip is assumed to be continuous everywhere 
inside the ellipse, together with its first derivatives. 

It is easy to see that the four products in IjBlOfl and 
their normal derivatives have the following continuity 
properties across the segment (—a, a): 



F 


F continuous? 


d n F continuous? 


ce n { 


<j>,q)Je n (z,q) 


Yes 


Yes 


se n ( 


(j), q)Jo n (z,q) 


Yes 


Yes 


ce„( 


(p,q)Ne n (z,q) 


Yes 


No 


se„( 


(j), q)No n (z,q) 


No 


Yes 



For a particle-in-the-box problem, therefore, the solu- 
tion must be of the form 



f = 2_j a nce n ((t), q)Je n (z, q) + ^ 7„se„(0, q)Jo n (z, q). 

n—0 n—1 

Coefficients a n and 7„ are determined by the boundary 
conditions imposed on the exterior boundary of the sys- 
tem. If the boundary of the ellipse is described by z = R, 
then the eigenstates can be labelled by two integer in- 
dices: the angular index n, the radial index j, and their 
parity under <f> — ► — <fi. More explicitly, the eigenfunctions 
are (up to normalization factors) 

^)(c(>,z)=ce n (<p,qe j )Je n (4>,q e nj ) , (Bll) 
^nj( r P, z ) = se n ((t),q° :j )Jo n ((t),q° :j ) , (B12) 

where q^j (or q°j) is defined as the j-th solution of 
Je n (R,q) = (or Jo n (R,q) = 0). The corresponding- 
energy eigenvalues are 

<? = !-fc 2 = !-%- • (B13) 
" J 2m 2m a 

2. Schrodinger particle in the field of two 
Aharonov-Bohm half-fluxes 

a. Boundary condition on a segment z = 

To solve the Schrodinger equation 

(— iV — eA) 2 ip + k 2 ip = (B14) 

with the vector potential A(r) describing two Aharonov- 
Bohm fluxes at positions (a,0) and (— a, 0), we perform 
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singular gauge transformation 



V> -> / exp ^ie y" A(r) • dr^j . (B15) 

Since the integral in the exponent depends on the path of 
integration, we need to introduce a branch cut in order 
to obtain single-valued functions. We chose the segment 
(—a, a) of the real axis as the branch cut: 



ip -> f exp 



(B16) 



where — ir < 9i < ir and < 62 < 2ir. This transfor- 
mation eliminates the vector potential, and the resulting 
equation is just the Helmholtz equation 

V 2 / + fc 2 / = , 

except that the solutions must have a branch cut on a 
segment (—a, a): 

f(x, y + e) = -f(x, y - e) for x £ (-a, a) . 

Elliptic coordinates provide a natural framework for im- 
posing the boundary condition: 

f(z = o,4>) = -f(z = o,-4>) , 

i.e. at z = the solution f(z = 0, <ft) must be an odd 
function of cf>. 

Starting from a general form 

/ = ^ ce n (4>)[a n Je n (z) + (3 n Ne n (z)\ 

n=0 

+ J2se n ((j))[ ln Jo n (z) + 5 n No n (z)] , (B17) 

71=1 

we obtain for f(z — 0, (f): 

(<j>)No n (0), 

n— n— 1 

which implies the following relation between a n and /3„: 




FIG. 12: Two half-fluxes inside an elliptical box. 



we obtain 



ln Jo' n {0) + S n No' n (0) = 



(B20) 



Equations (|B18I) and (|B20|) express the continuity of the 
wavefunction ip(x, y) and its normal derivative at the seg- 
ment (—a, a). As a result, the general expression for the 
wave function has the form 



n=0 



f( z , <t>) = X] an(q)ce n (<f>, q)Le n (z, q) 

^2'y n (q)se n ((f),q)Lo n (z,q) . (B21) 



To compactify the notation, we have introduced the fol- 
lowing functions: 



Le n (z,q) 



Je n (z,q) _ Ne n (z,q) 

Je n (0,q) Ne n (0,q) 

T , s _ Jo n {z,q) No n (z,q) 
Lo n (z,q)- jQ , Mq) 



(B22) 
(B23) 



At this point one can pose and solve several problems 
with various external boundary conditions, which we will 
consider next. 



a n Je n (0) + /3 n Ne n (0) = . (B18) 

Similarly, the requirement for ip(x,y) to have continuous 
normal derivative across the branch cut yields 

d d 

Q^f{Z,4>)\z=0 = faf(z,-<l>)\z=0 ■ 

Since j^f(z, (f>) at z — equals 
J2ce n (<t>)[anJe' n (0) + PnNe' n (0)] 

n=0 

+ J2^n(4>)hnJo' n (0) + S n No' n (0)] , (B19) 

n=l 



6. Schrodinger particle in an elliptical box in the presence 
of two half-fluxes. 

The wavefunctions are zero on the boundary of the 
elliptical box z = R, and therefore 

^ a n (q)ce n {4>, q)Le n {R, q) 
n=0 

+ ^2ln(q)se n {^,q)Lo n (R,q) = 0. (B24) 
n=l 

Using orthogonality of ce n ((j>) and se n (4>), we find that 
the eigenstates of are either even (labeled as "e" ) or odd 
(labeled as "o") functions of <f>, and their explicit form, 
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FIG. 13: Two half-fluxes: scattering problem. 



the the wavefunctions to the sum of the incident plane 
wave exp(ik • r) and the outgoing wave described by the 
second term. Together, they allow to determine the scat- 
tering amplitude unambiguously. To match the two ex- 
pressions, we write (|B30I in terms of Mathieu functions. 
Expansion of a plane wave e lk r in Mathieu functions ba- 
sis has the following form: 



^ik{x cos 6+y sin 6) 



^2 Pn(q)Ce n (z)ce n (cf))ce n {9) 



71=0 



cr n {q)Se n {z)se n {4>)se n {6) , (B31) 



where the coefficients p n (q) , a n (q) are related to the 
known standard factors p n (q) , s n (q) from the theory of 
Mathieu functions (cf. Table IV|l as 



up to normalization factors, is given by 



ce Tt 



where q^j and g° • are j-th solutions of 
Je n (R,q) Ne n (R,q) 



Je n (0,q) Ne n (0,q) 
JOn{R,q) _ No n {R,q) 
Jo' n {Q,q) No' n (0,q) 



, 




(B25) 
(B26) 

(B27) 
(B28) 



respectively. Therefore, just as in the case of the no- 
flux problem of the previous section, the eigenstates are 
classified according to parity under (j) — > —(f), and two 
integers (j, n), where j = 1, 2, 3, . . ., and n = (0), 1,2,... 
for even (odd) solutions. The energy of the eigenstates 
\n,j,o(e)) is 



E nj ~ 



H 2 4g 



o(e) 



2m 



2m 



(B29) 



c. Scattering of a Schrddinger particle by two fluxes 

Now we consider the problem of a particle scattered by 
two fluxes 4>o/2 and — 0o/2. To find the scattering cross- 
section er(n, n') and the scattering amplitude F(h, n') 
we need to match general solution f(z,<p;q) given by 
(|B21|) with the boundary conditions at infinity that con- 
tains only the incoming and the scattered waves: 

/n = e lkr + ^Pp^F(n,n') . (B30) 

Here, n = (cos 9, sin 9) specifies the direction of the in- 
cident wave and n' = (cos 0, sin 0) describes the direc- 
tion of scattered wave. The first condition, Eq. I|B21|I . 
encodes the information about the presence of two half- 
fluxes, while the second condition, Eq. (|B30|) . restricts 



p2r, 



0277 



P2r, 
1 

S2« 



&2n+l 



P2n+1 

i 

S2n+1 



(B32) 
(B33) 



Now we turn to the scattered wave: the most general 
expression for the solution, which contains only the out- 
going wave, is given by 



ikr 



F(n, n') 



Y / ^nCe n (^)ce n (9)He ( n 1 \z) 

{4,)se n {e)Ho^{z) . (B34) 



n=0 



71=1 



Functions He^z) and ffo'(z), which play the role of the 
Hankel functions in the theory of Mathieu equation, are 
defined as 



He\(z) = Je n (z)+iNe n (z) 
Ho^(z) = Jo n (z) + iNo n (z) 



Combining (|B2T|) . (fB3T|l and lB34l we find 

Je„(0, q) 



A n — — 



Pn 



Je n (0,q)+iNe n {0,q) 

Jo' n (Q,q) 

Jo' n (0,q) + iNo' n (0,q) t 



Pn 



(B35) 

(B36) 
(B37) 



We remind the reader that all quantities on the right 
hand side of (|B37|) are known, and to obtain the scatter- 
ing amplitude F(n, n') from l|B34() we only need to use 
the long distance expressions for He n and Ho n that are 
easy to find from Table IVl 



PnHdr. 



a n Ho„ 



9 

ir z q 



1/4 



At large distances, we have e z = 2r/a, the elliptic co- 
ordinate <f> is reduced to the ordinary polar angle, and 
therefore 



p n He n rj a n Ho n « — — 
\ir 2 q 



1/4 



— 27r/4 ^ ikr 



(B38) 
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Finally, substitution of these expressions into (|B34|) 
yields the following exact expression for the scattering 
amplitude: 



F(e, 



2 \ 1/4 



ir 2 q 



e" l,r/4 x 



£ 



Je n (0,q) 
He n {0,q) 



ce n ((f), q)ce n (9,q) 



E g ; ( ( °o' i g g ) ) se "^g) se "^g)) • ( B39 ) 



3. Dirac equation 

a. Relation between Schrodinger and Dirac problems 

Now we apply the technique we developed in this ap- 
pendix to the problem of a Dirac particle moving in the 
field of two half-fluxes, without specifying the external 
boundary condition at this point. The Hamiltonian Hq 
describing such a system is 

I { Px - eA x ) - i( Py - eA y )\ 

\(p x - eA x ) + i(p y - eA y ) J ' 

where (p x ,p y ) is the momentum operator and A is the 
vector potential of the two half-fluxes. Just as for the 
Schrodinger problem, we perform a unitary transforma- 
tion (|B16|) , which eliminates the vector potential from the 
Hamiltonian, but at the expense of introducing a branch 
cut. After choosing a segment of the real axis (—a, a) as 
the location of the branch cut, the eigenfunctions of the 
Hamiltonian 

H = e( ° Px ~ iPy ] 

\Px +iPy J 

can be found from the eigenfunctions of the Schrodinger 
problem, which we described in the previous section: 
note that if / = (u, v) T is an eigenfunction of H, then it 
is also an eigenfunction of 



pI+pI 



V Pi+PyJ 

Thus, both u and v are solutions of the wave equation 

(V 2 + E 2 )u = , (B40) 

and any solutions of our eigenvalue problem 

(px ~ ip y )v = Eu , (B41) 

(B42) 

, | , (B43) 

e(Px +ip y )u(x,y)l 




where u(r) is a solution of l|B40(l . The opposite statement 
is also clearly valid: for any u, which satisfies <|B40[1 . 
wavefunction 1B43|1 is an eigenfunction of Hamiltonian H 
with energy E. Thus, it naively appears that all solutions 
of the Dirac equation can be written as l|B43(l . We will 
return to the validity of this statement in a moment, but 
let us first apply (|B43|) to a simpler problem of a single- 
flux problem. 



b. The single-flux problem 

A problem of Schrodinger particle moving in the field 
of a single half-flux located at the origin is most con- 
veniently solved by eliminating the vector potential by 
means of ip{ v ) = ex P(*$/2)/(r) transformation, where 6 
is the polar angle. The resulting eigenvalue problem 



(V 2 +E 2 )f = 



(B44) 



must be solved in a space of wavefunctions /(r, 9) that 
have a branch-cut extending from to infinity. The cut 
can be chosen arbitrarily, e.g., as a straight line (0, +oo), 
or (—oo,0). The solutions are easily obtained in polar 
coordinates as 



/=£« 



,»(m+l/2)0 



(B45) 

Consider now a class of problems where the boundary 
conditions allow the particles to reach the origin - this 
excludes problems where impenetrable walls of finite ra- 
dius surround the flux. 

The requirement of square-integrability of the wave- 
functions requires that all coefficients d m , m > and 
c m ,m < must be set to zero except, possibly, do and 
C—i. The eigenfunctions corresponding to c_i and do are 
divergent at the origin, but square integrable. Ordinarily, 
physical realization of the flux requires the wavefunctions 
to be not only square-integrable, but also finite at r = 0: 
this eliminates the remaining arbitrariness and results in 
d m = for all m > 0, and c m = for all m < 0. Thus, for 
every angular channel m there is only one radial solution: 



- JO/2 f _ 



/=£< 



1 J\m-i/2\(Er) 



however, (|B44|I is used as an auxiliary tool to obtain solu- 
tions of the Dirac problem via Eq. I|B43(I . requirement of 
the wavefunctions being finite at r = is too restrictive: 
even if the upper components u of solutions are chosen 
to be finite, then the lower component of at least some 
wavefunctions will necessarily be divergent, but square 
integrable at the origin. To consider the most general 
situation, however, we choose u is the form of l|B45(l al- 
lowing at this stage all square integrable wavefunctions. 
Using l|B43() and the following expressions: 



d x ± id y = e ±i& | d r ± -dt 



(B46) 
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the solutions of the Dirac equation can be written as 



where 
Y (±) = 

Am 



±(m+i) 



J ±{m+\)( Er ) 



(m+i) 



(B47) 

Using the standard identities for Bessel functions, / can 
be written as 



a i(m+i)6> 



ie ie J m+i {Er)) 
1 \-ie^J_ m s(Er) 



(B48) 



For all values of m > (m < —2) the requirement of 
square-integrability demands that d m = (c m =0). At 
m = — 1, however there is an ambiguity. It is impossible 
a priori to decide which of the two radial functions 




^(Er) 
-ie ie J_i(Er) 



(B49) 



should be used. Either the upper or the lower compo- 
nent of the wavefunction is divergent, but still square 
integrable. The attempt to set both c_i and d-i to zero 
leads to the loss of completeness in the angular basis: the 
wavefunctions described by / cx exp(— i6/2) or, equiva- 
lently, the original gauge wavefunctions ip which do not 
have angular dependence, would be left out of the Hilbcrt 
space. On the other hand, to require that both solutions 
are present in the spectrum independently would be too 
much: the resulting set of basis functions is then over- 
complete. As was shown by Jackiw and Gerbert, the 
solution is to use one linear combination of the two solu- 
tions. Different regularizations of the problem then cor- 
respond to different choice of the linear combination. Im- 
portantly, not all linear combinations are mathematically 
allowed, and can emerge from the regularized problems: 
the relative phase of c_i and d-\ turns out to be fixed, 
and different boundary conditions, forming different self- 
adjoint extensions of the problem are described through 
a single parameter 9. For a given 9 the divergent part of 
the wavefunctions is 



sin u 

-ie 1 ^ cos ( 



(B50) 



Note that in this problem the basis functions have a use- 
ful property: both the upper component and the lower 
component can be simultaneously written as separable 
functions, i.e. a product of two functions that depend 
only on 9 and only on r. This is a peculiarity of the 
single-flux problem, and in general this property does not 



hold. Nevertheless, although the analysis will be slightly 
more complicated, all essential properties of the single- 
flux problem such as the necessity of additional boundary 
conditions at flux locations and the form of the boundary 
conditions remain valid in the case of two-flux problem 
as well. 



Dirac particle in the presence of two half-fluxes 



Once again, we use (|B43|) to construct solutions of 
Dirac equation from the solutions (|B21Jl of the wave equa- 
tion in the presence of two half-fluxes. In elliptical coor- 
dinates, 



d x + id v 



i sinh z cos < 



i cosh z sin < 



The eigenfunctions of the Dirac equation obtained from 
(|B43|) and l|B21|l can be written now as 



n=0 



M<z)xl +) + y>«(<7)xl _) 



n=l 



(B51) 



where spinors x!^ are equal to 



A+) 



CG n (c/>)Le n (z,q) 



{4>)Le' n {z) 
E a sinh z cos c 



■ice' n (<fi)Le n (z) 
- a cosh z sin <p 



(B52) 



and 



A-) 



se n ((j))Lo n (z,q) 



j se n (4>)Lo' n (z) 
E a sinh z cos < 



i se' n (4>)Lo n (z) 
- a cosh z sin (f> 



(B53) 



So far we almost literally followed the route of the single 
half-flux problem. However, even a cursory examination 

of the solutions Xn shows that something is amiss. In 
obtaining l|B21(l and then Xn \ we never discarded any 

solutions, and yet the upper components of all Xn^ are 
always perfectly regular. So, where are the wavefunctions 
with upper components divergent near the half-fluxes? 

Before we answer this question, let us examine the 
lower components of the solutions. The elliptical coor- 
dinates of the flux located at (x,y) = (a,0) are (z,<j>) — 
(0, 0). Therefore, in the vicinity of this flux, we have 



1 



1 



a sinh z cos (f> — a cosh z sin <f> a(z — icf>) ^2ap\ ' 

(B54) 

where p\ is the distance between the point (x, y) and 
the flux, and 9\ G (— 7r, 7r) is the polar angle shown in 
Fig- El Similarly, near the second flux at (— a,0), we 
have (z, (f) w (0,7r), and therefore 



1 



1 



a sinh z cos < 



acoshzsin</> a(z — iq 



\J2api 
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Since Le n (z = 0, q) =0 and Lo' n (z = 0, q) = 0, in the 
vicinity of the first flux, the divergent part of Xn^ is 



X 



X 



(+) = ice n (0,q)Le' n (0 : q) e^/ 2 / o" 
E ^2b~p~l \ l y 

( _j _ se' n (0,q)Lo n (0,q) e^' 2 ( (A 



E 



V2api \ 1 



(B55) 
(B56) 



The divergent part of the wavefunctions near the sec- 
ond flux is almost identical: 

(+) _ ice B (7r, g )L<(0, g ) e ^/='-W2 / \ 

Xn ~ e 7^r[i) ' ( j 
x - 1 T^rW • (B58) 

Comparison with the single half-flux boundary condition 
at the flux location (|B50I) suggests that the wavefunctions 
Xn, which we just constructed, are merely one of many 
possible self-adjoint extensions. This extension contains 
wavefunctions which near both fluxes have regular upper 
components and divergent lower component. 

What about the other self-adjoint extensions? Where 
did we lose them? After all, we only repeated the steps 
for the single flux Aharonov-Bohm problem, where this 
approach allowed us to find all self-adjoint extensions. 
Why didn't they naturally appear from l|B43j) ? 

To understand what went wrong, consider again solu- 
tions of the Schrodinger equations for the single half- 
flux (|B45|) and compare them to general solution of 
the Schrodinger solution for the two half-fluxes problem 
(IB21|I . One glaring distinction is that the former con- 
tains solutions divergent near the fluxes, while the latter 
does not! To be sure, there are certainly solutions of the 
two-flux problem that diverge near the fluxes - more- 
over, we constructed them explicitly, when we found the 
lower components of Xn ■ The paradox appeared because 
these divergent solutions do not have a separable form 
/i (z)/2(0) because the fluxes are point-like defects. The 
line connecting the two fluxes, and the elliptical coordi- 
nates used to describe it were just tools, which allowed 
us to analyze the problem analytically, but the under- 
lying structure of the divergencies is still that of point 



like defects. They are, naturally, most easily described 
in local polar coordinate system around each flux. The 
functional form of the singularities, which has a struc- 
ture Fi(pi)F2(0i), is incompatible with a factorization 
fi( z )f2((f>), as can be seen for example from (JB54I1 . 

In short, we were able to find the self-adjoint exten- 
sion Eq. (|B51|) - (|B53|) analytically in elliptical coordi- 
nates so easily only because we started with the regular 
upper components. Then it was easy to construct a com- 
plete basis of wavefunctions of the form ce n ((j))Le n (z) and 
se n (4>)Lo n (z). Only then we worked out the lower com- 
ponents, which ended up divergent near flux locations, 
but not factorizable. 

Can we find other self-adjoint extensions using ellipti- 
cal coordinates? There is one other case that is easy to 
solve: wavefunctions with the regular lower component. 
Rather than using l|B43j) . we could have tried the form 



/ = 



^(Px - ip y )v(x,y) 
v{x,y) 



(B59) 



We would start from the regular solutions of the wave 
equation for v. 

v(<j>, z) = a n ce n ((/))Le n (z) + (3 n se n (<j))Lo n (z) , 

and then would apply operator (j> x — ip y )/E to find the 
upper components of the spinors forming the basis. This 
naturally would result in a self-adjoint extension with 
the regular lower components, and divergent upper com- 
ponents - the opposite of the first self-adjoint extension 
we found in l|B51IB53(l . Note that these two self-adjoint 
extensions just found, with the regular upper or lower 
components, are the ones that follow from the simplest 
physical regularizations^ of the problem obtained by re- 
placing each of the the infinitely thin Aharonov-Bohm 
strings by a solenoid of finite radius. 

Can one do better and construct the basis analytically 
for the general case, with arbitrary parameters #1,2 char- 
acterizing via l|B50|) each of the fluxes? At present, we 
do not know the answer to this question, and leave this 
interesting problem for future study. 
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